Hello All,
Not sure if this is stepping on toes, but here is what I thought of to deal
with pivoting. I would only need to modify the constructor of the
QRDecompImpl class:
public QRDecompositionPivotImpl(RealMatrix matrix) {
final int m = matrix.getRowDimension();
final int n = matrix.getColumnDimension();
final int mn = FastMath.min(m, n);
qrt = matrix.transpose().getData();
rDiag = new double[FastMath.min(m, n)];
cachedQ = null;
cachedQT = null;
cachedR = null;
cachedH = null;
/*
* The QR decomposition of a matrix A is calculated using
Householder
* reflectors by repeating the following operations to each minor
* A(minor,minor) of A:
*/
pivots = new int[n];
for (int i = 0; i < n; i++) {
pivots[i] = i;
}
int pivotIdx = -1;
double bestNorm = 0.0;
for (int minor = 0; minor < mn; minor++) {
bestNorm = 0.0;
pivotIdx = 0;
for (int i = minor; i < n; i++) {
final double[] qrtMinor = qrt[i];
double xNormSqr = 0.0;
for (int row = minor; row < m; row++) {
final double c = qrtMinor[row];
xNormSqr += c * c;
}
if (xNormSqr > bestNorm) {
bestNorm = xNormSqr;
pivotIdx = i;
}
}
if( pivotIdx != minor){
int pvt = pivots[minor];
pivots[minor] = pivots[pivotIdx];
double[] qswp = qrt[minor];
qrt[minor] = qrt[pivotIdx];
for( int i = minor + 1; i < n; i++){
if( i <= pivotIdx ){
final int itmp = pivots[i];
pivots[i] = pvt;
pvt = itmp;
final double[] tmp = qrt[i];
qrt[i] = qswp;
qswp=tmp;
}
}
}
final double[] qrtMinor = qrt[minor];
/*
* Let x be the first column of the minor, and a^2 = |x|^2.
* x will be in the positions qr[minor][minor] through
qr[m][minor].
* The first column of the transformed minor will be (a,0,0,..)'
* The sign of a is chosen to be opposite to the sign of the
first
* component of x. Let's find a:
*/
final double a = (qrtMinor[minor] > 0) ?
-FastMath.sqrt(bestNorm) : FastMath.sqrt(bestNorm);
rDiag[minor] = a;
if (a != 0.0) {
/*
* Calculate the normalized reflection vector v and
transform
* the first column. We know the norm of v beforehand: v =
x-ae
* so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
* a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
* Here <x, e> is now qr[minor][minor].
* v = x-ae is stored in the column at qr:
*/
qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
/*
* Transform the rest of the columns of the minor:
* They will be transformed by the matrix H = I-2vv'/|v|^2.
* If x is a column vector of the minor, then
* Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2
v.
* Therefore the transformation is easily calculated by
* subtracting the column vector (2<x,v>/|v|^2)v from x.
*
* Let 2<x,v>/|v|^2 = alpha. From above we have
* |v|^2 = -2a*(qr[minor][minor]), so
* alpha = -<x,v>/(a*qr[minor][minor])
*/
for (int col = minor + 1; col < n; col++) {
final double[] qrtCol = qrt[col];
double alpha = 0;
for (int row = minor; row < m; row++) {
alpha -= qrtCol[row] * qrtMinor[row];
}
alpha /= a * qrtMinor[minor];
// Subtract the column vector alpha*v from x.
for (int row = minor; row < m; row++) {
qrtCol[row] -= alpha * qrtMinor[row];
}
}
}
}
}
All I am doing is looking forward and taking the next column with the
largest norm. Then I rearrange the Q's. Is this what you had in mind Chris?
The result will be Q,R and the pivot array for which we can implement a
getter?
-Greg
