mikera.matrixx.decompose.impl.hessenberg.TridiagonalDecompositionHouseholder Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of vectorz Show documentation
Show all versions of vectorz Show documentation
Fast double-precision vector and matrix maths library for Java, supporting N-dimensional numeric arrays.
/*
* Copyright (c) 2009-2013, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package mikera.matrixx.decompose.impl.hessenberg;
import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.decompose.impl.qr.QRHelperFunctions;
/**
*
* Performs a TridiagonalSimilarDecomposition similar tridiagonal decomposition on a square symmetric input matrix.
* Householder vectors perform the similar operation and the symmetry is taken advantage
* of for good performance.
*
*
* Finds the decomposition of a matrix in the form of:
*
* A = O*T*OT
*
* where A is a symmetric m by m matrix, O is an orthogonal matrix, and T is a tridiagonal matrix.
*
*
* This implementation is based off of the algorithm described in:
*
* David S. Watkins, "Fundamentals of Matrix Computations," Second Edition. Page 349-355
*
*
* @author Peter Abeles
*/
public class TridiagonalDecompositionHouseholder {
/**
* Only the upper right triangle is used. The Tridiagonal portion stores
* the tridiagonal matrix. The rows store householder vectors.
*/
private Matrix QT;
// The size of the matrix
private int N;
// temporary storage
private double w[];
// gammas for the householder operations
private double gammas[];
// temporary storage
private double b[];
public TridiagonalDecompositionHouseholder() {
N = 1;
w = new double[N];
b = new double[N];
gammas = new double[N];
}
/**
* Returns the internal matrix where the decomposed results are stored.
* @return
*/
public AMatrix getQT() {
return QT;
}
public void getDiagonal(double[] diag, double[] off) {
for( int i = 0; i < N; i++ ) {
diag[i] = QT.data[i*N+i];
if( i+1 < N ) {
off[i] = QT.data[i*N+i+1];
}
}
}
/**
* Extracts the tridiagonal matrix found in the decomposition.
*
* @return The extracted T matrix.
*/
public AMatrix getT() {
Matrix T = Matrix.create(N,N);
T.data[0] = QT.data[0];
for( int i = 1; i < N; i++ ) {
T.set(i,i, QT.get(i,i));
double a = QT.get(i-1,i);
T.set(i-1,i,a);
T.set(i,i-1,a);
}
if( N > 1 ) {
T.data[(N-1)*N+N-1] = QT.data[(N-1)*N+N-1];
T.data[(N-1)*N+N-2] = QT.data[(N-2)*N+N-1];
}
return T;
}
/**
* An orthogonal matrix that has the following property: T = QTAQ
*
* @return The extracted Q matrix.
*/
public AMatrix getQ( boolean transposed ) {
Matrix Q = Matrix.createIdentity(N);
for( int i = 0; i < N; i++ ) w[i] = 0;
if( transposed ) {
for( int j = N-2; j >= 0; j-- ) {
w[j+1] = 1;
for( int i = j+2; i < N; i++ ) {
w[i] = QT.data[j*N+i];
}
QRHelperFunctions.rank1UpdateMultL(Q,w,gammas[j+1],j+1,j+1,N);
}
} else {
for( int j = N-2; j >= 0; j-- ) {
w[j+1] = 1;
for( int i = j+2; i < N; i++ ) {
w[i] = QT.get(j,i);
}
QRHelperFunctions.rank1UpdateMultR(Q,w,gammas[j+1],j+1,j+1,N,b);
}
}
return Q;
}
/**
* Decomposes the provided symmetric matrix.
*
* @param A Symmetric matrix that is going to be decomposed. Not modified.
*/
public boolean decompose( AMatrix A ) {
init(A);
for( int k = 1; k < N; k++ ) {
similarTransform(k);
}
return true;
}
/**
* Computes and performs the similar a transform for submatrix k.
*/
private void similarTransform( int k) {
double t[] = QT.data;
// find the largest value in this column
// this is used to normalize the column and mitigate overflow/underflow
double max = 0;
int rowU = (k-1)*N;
for( int i = k; i < N; i++ ) {
double val = Math.abs(t[rowU+i]);
if( val > max )
max = val;
}
if( max > 0 ) {
// -------- set up the reflector Q_k
double tau = QRHelperFunctions.computeTauAndDivide(k,N,t,rowU,max);
// write the reflector into the lower left column of the matrix
double nu = t[rowU+k] + tau;
QRHelperFunctions.divideElements(k+1,N,t,rowU,nu);
t[rowU+k] = 1.0;
double gamma = nu/tau;
gammas[k] = gamma;
// ---------- Specialized householder that takes advantage of the symmetry
householderSymmetric(k,gamma);
// since the first element in the householder vector is known to be 1
// store the full upper hessenberg
t[rowU+k] = -tau*max;
} else {
gammas[k] = 0;
}
}
/**
* Performs the householder operations on left and right and side of the matrix. QTAQ
* @param row Specifies the submatrix.
*
* @param gamma The gamma for the householder operation
*/
public void householderSymmetric( int row , double gamma )
{
int startU = (row-1)*N;
// compute v = -gamma*A*u
for( int i = row; i < N; i++ ) {
double total = 0;
// the lower triangle is not written to so it needs to traverse upwards
// to get the information. Reduces the number of matrix writes need
// improving large matrix performance
for( int j = row; j < i; j++ ) {
total += QT.data[j*N+i]*QT.data[startU+j];
}
for( int j = i; j < N; j++ ) {
total += QT.data[i*N+j]*QT.data[startU+j];
}
w[i] = -gamma*total;
}
// alpha = -0.5*gamma*u^T*v
double alpha = 0;
for( int i = row; i < N; i++ ) {
alpha += QT.data[startU+i]*w[i];
}
alpha *= -0.5*gamma;
// w = v + alpha*u
for( int i = row; i < N; i++ ) {
w[i] += alpha*QT.data[startU+i];
}
// A = A + w*u^T + u*w^T
for( int i = row; i < N; i++ ) {
double ww = w[i];
double uu = QT.data[startU+i];
int rowA = i*N;
for( int j = i; j < N; j++ ) {
// only write to the upper portion of the matrix
// this reduces the number of cache misses
QT.data[rowA+j] += ww*QT.data[startU+j] + w[j]*uu;
}
}
}
/**
* If needed declares and sets up internal data structures.
*
* @param A Matrix being decomposed.
*/
public void init( AMatrix A ) {
if( A.rowCount() != A.columnCount())
throw new IllegalArgumentException("Must be square");
if( A.columnCount() != N ) {
N = A.columnCount();
if( w.length < N ) {
w = new double[ N ];
gammas = new double[N];
b = new double[N];
}
}
QT = Matrix.create(A);
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy