TestHouseholderColQR.java

/*
 * Copyright (c) 2009-2013, Peter Abeles. All Rights Reserved.
 *
 * This file is part of Efficient Java Matrix Library (EJML).
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package mikera.matrixx.decompose.impl.qr;

import java.util.Random;

import mikera.matrixx.Matrix;
import mikera.matrixx.algo.Multiplications;
import mikera.matrixx.impl.AStridedMatrix;

import org.junit.Test;

import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertTrue;

/**
 * @author Peter Abeles
 */
public class TestHouseholderColQR extends GenericQrCheck {

    Random rand = new Random(0xff);
    
    @Override
    protected QRDecomposition createQRDecomposition(boolean compact)
    {
        return new HouseholderColQR(compact);
    }

    /**
     * Internal several householder operations are performed.  This
     * checks to see if the householder operations and the expected result for all the
     * submatrices.
     */
    @Test
    public void householder() {
        int width = 5;

        for( int i = 0; i < width; i++ ) {
            checkSubHouse(i , width);
        }
    }

    private void checkSubHouse(int w , int width) {
        DebugQR qr = new DebugQR(width,width);

        Matrix A = Matrix.createRandom(width, width);

        qr.householder(w,A);

//        SimpleMatrix U = new SimpleMatrix(width,1, true, qr.getQR()[w]).extractMatrix(w,width,0,1);
        Matrix temp = Matrix.create(width,1);
        temp.setElements(qr.getQR()[w]);
        AStridedMatrix U = temp.subMatrix(w, width-w, 0, 1);
        
        U.set(0,0,1); // this is not explicity set and is assumed to be 1
        Matrix I = Matrix.createIdentity(width-w);
//        SimpleMatrix Q = I.minus(U.mult(U.transpose()).scale(qr.getGamma()));
        Matrix temp1 = Multiplications.multiply(U, U.getTranspose());
        temp1.scale(qr.getGamma());
        I.sub(temp1);
        Matrix Q = I;


        // check the expected properties of Q
        assertTrue(Q.epsilonEquals(Q.getTranspose(),1e-6));
        assertTrue(Q.epsilonEquals(Q.inverse(),1e-6));

//        SimpleMatrix result = Q.mult(A.extractMatrix(w,width,w,width));
        AStridedMatrix result = Multiplications.multiply(Q, A.subMatrix(w,width-w,w,width-w));

        for( int i = 1; i < width-w; i++ ) {
            assertEquals(0,result.get(i,0),1e-5);
        }
    }

    /**
     * Check the results of this function against basic matrix operations
     * which are equivalent.
     */
    @Test
    public void updateA() {
        int width = 5;

        for( int i = 0; i < width; i++ )
            checkSubMatrix(width,i);
    }

    private void checkSubMatrix(int width , int w ) {
        DebugQR qr = new DebugQR(width,width);

        double gamma = 0.2;
        double tau = 0.75;

        Matrix U = Matrix.createRandom(width, 1);
        Matrix A = Matrix.createRandom(width, width);

        qr.convertToColumnMajor(A);

        // compute the results using standard matrix operations
        Matrix I = Matrix.createIdentity(width-w);

//        SimpleMatrix u_sub = U.extractMatrix(w,width,0,1);
        AStridedMatrix u_sub = U.subMatrix(w, width-w, 0, 1);
        u_sub.set(0,0,1);// assumed to be 1 in the algorithm
//        SimpleMatrix A_sub = A.extractMatrix(w,width,w,width);
        AStridedMatrix A_sub = A.subMatrix(w,width-w,w,width-w);
//        SimpleMatrix expected = I.minus(u_sub.mult(u_sub.transpose()).scale(gamma)).mult(A_sub);
        Matrix temp1 = Multiplications.multiply(u_sub, u_sub.getTranspose());
        temp1.scale(gamma);
        I.sub(temp1);
        Matrix expected = Multiplications.multiply(I, A_sub);

        qr.updateA(w,U.asDoubleArray(),gamma,tau);

        double[][] found = qr.getQR();

        for( int i = w+1; i < width; i++ ) {
            assertEquals(U.get(i,0),found[w][i],1e-8);
        }

        // the right should be the same
        for( int i = w; i < width; i++ ) {
            for( int j = w+1; j < width; j++ ) {
                double a = expected.get(i-w,j-w);
                double b = found[j][i];

                assertEquals(a,b,1e-6);
            }
        }
    }

    private static class DebugQR extends HouseholderColQR
    {

        public DebugQR( int numRows , int numCols ) {
            super(false);
            error = false;

            this.numCols = numCols;
            this.numRows = numRows;
            minLength = Math.min(numRows,numCols);
            int maxLength = Math.max(numRows,numCols);

            dataQR = new double[ numCols ][  numRows ];
            v = new double[ maxLength ];
            gammas = new double[ minLength ];
            
            this.numCols = numCols;
            this.numRows = numRows;
        }

        public void householder( int j , Matrix A ) {
            convertToColumnMajor(A);

            super.householder(j);
        }

        protected void convertToColumnMajor(Matrix A) {
            super.convertToColumnMajor(A);
        }

        public void updateA( int w , double u[] , double gamma , double tau ) {
            System.arraycopy(u,0,this.dataQR[w],0,u.length);
            this.gamma = gamma;
            this.tau = tau;

            super.updateA(w);
        }

        public double getGamma() {
            return gamma;
        }
    }
}