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 Sun, 07 Aug 2011 16:45:33 GMT
Your selection methodology might be more complicated than pairwise
comparison. So I think you would have to pass the entire set of remaining
columns. That is an implementation detail, however.

On the interface, yes the public face should not expose the constructor
which takes the PivotStrategyComparator. It should be part of a protected
constructor or, I suppose, we could have a static factory class which spits
out different implementations.

-Greg

On Sun, Aug 7, 2011 at 11:33 AM, Chris Nix <chris.nix@gmail.com> wrote:

> 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