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
>>
>
>
