commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Chris Nix <chris....@gmail.com>
Subject Re: [math] Pivoting in QR decomposition
Date Sun, 07 Aug 2011 12:28:11 GMT
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 <gsterijevski@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
>

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