commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Greg Sterijevski <gsterijev...@gmail.com>
Subject Re: [math] Pivoting in QR decomposition
Date Sat, 06 Aug 2011 20:34:12 GMT
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

Mime
  • Unnamed multipart/alternative (inline, None, 0 bytes)
View raw message