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....@gmail.com> wrote:
> Oooops, the ********* below should read MATH-630.
>
> Chris
>
> On 7 August 2011 13:28, Chris Nix <chris....@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 <gsterijev...@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
> >>
> >
> >
>
