mikera.matrixx.decompose.impl.bidiagonal.BidiagonalRow 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.bidiagonal;
import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.decompose.IBidiagonalResult;
import mikera.matrixx.decompose.impl.qr.QRHelperFunctions;
/**
*
* Performs a {@link org.ejml.alg.dense.decomposition.bidiagonal.BidiagonalDecomposition} using
* householder reflectors. This is efficient on wide or square matrices.
*
*
* @author Peter Abeles
*/
public class BidiagonalRow {
// A combined matrix that stores te upper Hessenberg matrix and the orthogonal matrix.
private Matrix UBV;
// number of rows
private int m;
// number of columns
private int n;
// the smaller of m or n
private int min;
// the first element in the orthogonal vectors
private double gammasU[];
private double gammasV[];
// temporary storage
private double b[];
private double u[];
private boolean compact;
private BidiagonalRow() {
}
/**
* Computes the decomposition of the provided matrix.
*
* @param A The matrix that is being decomposed. Not modified.
* @return an IBidiagonalResult object
*/
public static IBidiagonalResult decompose(AMatrix A) {
BidiagonalRow temp = new BidiagonalRow();
return temp._decompose(A, false);
}
/**
* Computes the decomposition of the provided matrix.
*
* @param A The matrix that is being decomposed. Not modified.
* @param compact if true, result matrices have zero-filled regions trimmed off.
* @return an IBidiagonalResult object
*/
public static IBidiagonalResult decompose(AMatrix A, boolean compact) {
BidiagonalRow temp = new BidiagonalRow();
return temp._decompose(A, compact);
}
/**
* Computes the decomposition of the provided matrix.
*
* @param A The matrix that is being decomposed. Not modified.
* @param compact If true, result matrices have zero-filled regions trimmed off
* @return If it detects any errors or not.
*/
private IBidiagonalResult _decompose(AMatrix A, boolean compact)
{
this.compact = compact;
UBV = Matrix.create(A);
m = UBV.rowCount();
n = UBV.columnCount();
min = Math.min(m, n);
int max = Math.max(m, n);
b = new double[max+1];
u = new double[max+1];
gammasU = new double[m];
gammasV = new double[n];
for( int k = 0; k < min; k++ ) {
// UBV.print();
computeU(k);
// System.out.println("--- after U");
// UBV.print();
computeV(k);
// System.out.println("--- after V");
// UBV.print();
}
return new BidiagonalRowResult(getU(), getB(), getV());
}
/**
* Returns the bidiagonal matrix.
*
* @return The bidiagonal matrix.
*/
private AMatrix getB() {
Matrix B = handleB(m,n,min);
//System.arraycopy(UBV.data, 0, B.data, 0, UBV.getNumElements());
B.set(0,0,UBV.get(0,0));
for( int i = 1; i < min; i++ ) {
B.set(i,i, UBV.get(i,i));
B.set(i-1,i, UBV.get(i-1,i));
}
if( n > m )
B.set(min-1,min,UBV.get(min-1,min));
return B;
}
private Matrix handleB( int m , int n , int min ) {
int w = n > m ? min + 1 : min;
if( compact ) {
return Matrix.create(min,w);
} else {
return Matrix.create(m,n);
}
}
/**
* Returns the orthogonal U matrix.
*
* @return The extracted Q matrix.
*/
private AMatrix getU() {
Matrix U = handleU(m,n,min);
for( int i = 0; i < m; i++ ) u[i] = 0;
for( int j = min-1; j >= 0; j-- ) {
u[j] = 1;
for( int i = j+1; i < m; i++ ) {
u[i] = UBV.get(i,j);
}
QRHelperFunctions.rank1UpdateMultR(U,u,gammasU[j],j,j,m,this.b);
}
return U;
}
private Matrix handleU( int m, int n , int min ) {
if( compact ){
return Matrix.createIdentity(m, min);
} else {
return Matrix.createIdentity(m, m);
}
}
/**
* Returns the orthogonal V matrix.
*
* @return The extracted Q matrix.
*/
private AMatrix getV() {
Matrix V = handleV(m,n,min);
// UBV.print();
// todo the very first multiplication can be avoided by setting to the rank1update output
for( int j = min-1; j >= 0; j-- ) {
u[j+1] = 1;
for( int i = j+2; i < n; i++ ) {
u[i] = UBV.get(j,i);
}
QRHelperFunctions.rank1UpdateMultR(V,u,gammasV[j],j+1,j+1,n,this.b);
}
return V;
}
private Matrix handleV( int m , int n , int min ) {
int w = n > m ? min + 1 : min;
if( compact ) {
return Matrix.createIdentity(n, w);
} else {
return Matrix.createIdentity(n, n);
}
}
protected void computeU( int k) {
double b[] = UBV.data;
// find the largest value in this column
// this is used to normalize the column and mitigate overflow/underflow
double max = 0;
for( int i = k; i < m; i++ ) {
// copy the householder vector to vector outside of the matrix to reduce caching issues
// big improvement on larger matrices and a relatively small performance hit on small matrices.
double val = u[i] = b[i*n+k];
val = Math.abs(val);
if( val > max )
max = val;
}
if( max > 0 ) {
// -------- set up the reflector Q_k
double tau = QRHelperFunctions.computeTauAndDivide(k,m,u ,max);
// write the reflector into the lower left column of the matrix
// while dividing u by nu
double nu = u[k] + tau;
QRHelperFunctions.divideElements_Bcol(k+1,m,n,u,b,k,nu);
u[k] = 1.0;
double gamma = nu/tau;
gammasU[k] = gamma;
// ---------- multiply on the left by Q_k
QRHelperFunctions.rank1UpdateMultR(UBV,u,gamma,k+1,k,m,this.b);
b[k*n+k] = -tau*max;
} else {
gammasU[k] = 0;
}
}
protected void computeV(int k) {
double b[] = UBV.data;
int row = k*n;
// find the largest value in this column
// this is used to normalize the column and mitigate overflow/underflow
double max = QRHelperFunctions.findMax(b,row+k+1,n-k-1);
if( max > 0 ) {
// -------- set up the reflector Q_k
double tau = QRHelperFunctions.computeTauAndDivide(k+1,n,b,row,max);
// write the reflector into the lower left column of the matrix
double nu = b[row+k+1] + tau;
QRHelperFunctions.divideElements_Brow(k+2,n,u,b,row,nu);
u[k+1] = 1.0;
double gamma = nu/tau;
gammasV[k] = gamma;
// writing to u could be avoided by working directly with b.
// requires writing a custom rank1Update function
// ---------- multiply on the left by Q_k
QRHelperFunctions.rank1UpdateMultL(UBV,u,gamma,k+1,k+1,n);
b[row+k+1] = -tau*max;
} else {
gammasV[k] = 0;
}
}
}© 2015 - 2025 Weber Informatics LLC | Privacy Policy