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:29:08 GMT
Oooops, the ********* below should read MATH-630.

Chris

On 7 August 2011 13:28, Chris Nix <chris.nix@gmail.com> wrote:

> 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