Thanks, Greg, for looking more at this. Apologies I've not been able to
focus on this too much recently.
Yes, the approach above works, however the solvers also require a change so
that they don't unexpectedly solve for a permuted matrix. Additionally, a
getRank() method can now be implemented much more efficiently than the
getRank() from SingularValueDecomposition.
I submitted an initial patch at *********** with these bits working, however
this patch introduces a new RRQRDecomposition interface extending
QRDecomposition. In light of the insights above from the team, I'll
implement it instead within the existing class structure and re-submit.
If the pivoting code is to be included in QRDecomposition, then perhaps an
extra constructor is required to perform column pivoting since the RRQR
decomposition of matrix A produces Q, R and P such A = QRP^T and would break
any existing code that expects a plain decomposition of A=QR.
Chris.
On 6 August 2011 21:34, Greg Sterijevski <gsterijev...@gmail.com> wrote:
> 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
>
