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 16:33:17 GMT
To avoid having the strategy method being overriden in the constructor, we
could instead implement a constructor that takes a column Comparator that
determines if two columns should be exchanged at each stage.

In the interest of maintaining a clean interface, I don't know if this
constructor should be public though.

Chris

On 7 August 2011 16:43, Greg Sterijevski <gsterijevski@gmail.com> wrote:

> Chris,
>
> In regard to the pivoting, do you think that it would be useful to allow
> subclasses to pivot using other strategies? The pivoting I have looks for
> the next column with the largest norm. For most garden variety problems
> this
> will be okay. However, you can (I am not sure how effective this will be)
> choose the column which has the lowest [average] inner product (ie is least
> correlated in some sense to the remaining columns).
>
> The easiest way to accomplish this is to present the remainder columns to
> some method, and have that method return an index... The method would be
> protected so that it could be overridden  by classes wishing to modify that
> behavior. The only problem I see is that the actual decomposition would
> need
> to be moved out of the constructor (we would have an overridable method
> being called in the constructor).
>
> -Greg
>
> On Sun, Aug 7, 2011 at 7:29 AM, Chris Nix <chris.nix@gmail.com> wrote:
>
> > 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