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 resubmit.
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 =
> xae
> * so v^2 = <xae,xae> = <x,x>2a<x,e>+a^2<e,e>
=
> * a^2+a^22a<x,e> = 2a*(a  <x,e>).
> * Here <x, e> is now qr[minor][minor].
> * v = xae 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 = I2vv'/v^2.
> * If x is a column vector of the minor, then
> * Hx = (I2vv'/v^2)x = x2vv'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
>
