source: public/Misc/HouseholderReflection.java @ 207

/*
 * Copyright (c) Numerical Method Inc.
 * http://www.numericalmethod.com/
 * 
 * THIS SOFTWARE IS LICENSED, NOT SOLD.
 * 
 * YOU MAY USE THIS SOFTWARE ONLY AS DESCRIBED IN THE LICENSE.
 * IF YOU ARE NOT AWARE OF AND/OR DO NOT AGREE TO THE TERMS OF THE LICENSE,
 * DO NOT USE THIS SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITH NO WARRANTY WHATSOEVER,
 * EITHER EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
 * ANY WARRANTIES OF ACCURACY, ACCESSIBILITY, COMPLETENESS,
 * FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABILITY, NON-INFRINGEMENT,
 * TITLE AND USEFULNESS.
 * 
 * IN NO EVENT AND UNDER NO LEGAL THEORY,
 * WHETHER IN ACTION, CONTRACT, NEGLIGENCE, TORT, OR OTHERWISE,
 * SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
 * ANY CLAIMS, DAMAGES OR OTHER LIABILITIES,
 * ARISING AS A RESULT OF USING OR OTHER DEALINGS IN THE SOFTWARE.
 */
package com.numericalmethod.suanshu.matrix.doubles.dense.operation.factorization;

import com.numericalmethod.suanshu.vector.doubles.Vector;
import com.numericalmethod.suanshu.matrix.doubles.Matrix;
import com.numericalmethod.suanshu.matrix.doubles.dense.DenseMatrix;
import com.numericalmethod.suanshu.matrix.doubles.dense.PermutationMatrix;
import com.numericalmethod.suanshu.matrix.doubles.dense.UpperTriangularMatrix;
import com.numericalmethod.suanshu.matrix.doubles.dense.operation.CreateMatrix;
import com.numericalmethod.suanshu.matrix.doubles.dense.operation.Householder;
import com.numericalmethod.suanshu.matrix.doubles.dense.operation.Householder.Context;
import com.numericalmethod.suanshu.misc.SuanShuUtils;
import com.numericalmethod.suanshu.vector.doubles.dense.DenseVector;
import static com.numericalmethod.suanshu.number.DoubleUtils.*;
import static com.numericalmethod.suanshu.vector.doubles.dense.operation.CreateVector.*;

/**
 * Successive Householder reflections gradually transform a matrix <code>A</code> to the upper triangular form.
 *
 * <p>
 * For example, the first step is to multiply <code>A</code> with the Householder matrix <code>Q<sub>1</sub></code> so that
 * matrix <code>Q1 %*% A</code> has zeros in the left column (except for the first row).
 * That is,
 *
 * <blockquote><code><pre>
 *             | a1 * ... * |
 *             | 0          |
 * Q1 %*% A =  | .          |
 *             | .     A'   |
 *             | .          |
 *             | 0          |
 * </pre></code></blockquote>
 *
 * At the end,
 * <blockquote><code><pre>
 * Qn %*% ... %*% Q1 %*% A = R
 * </pre></code></blockquote>
 * where <code>R</code> is upper triangular.
 *
 * <p>
 * Householder reflection has a better numerical stability than the Gram-Schmidt process.
 *
 * @author Haksun Li
 *
 * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections">Wikipedia: Using Householder reflections</a>
 */
public class HouseholderReflection implements QRDecomposition {

    /**
     * number of rows
     */
    public final int nrows;
    /**
     * number of columns
     */
    public final int ncols;
    /**
     * a precision parameter: when a number |x| ≀ ε, it is considered 0
     */
    public final double epsilon;
    private UpperTriangularMatrix R;//nCols x nCols
    private Vector[] cols;
    private Householder[] Hs;

    /**
     * Construct an instance of the Householder reflection process to orthogonalize a matrix.
     * 
     * @param A a matrix, where #rows >= #columns
     * @param epsilon a precision parameter: when a number |x| ≀ ε, it is considered 0
     * @throws IllegalArgumentException if #rows < #columns
     */
    public HouseholderReflection(Matrix A, double epsilon) {
        SuanShuUtils.assertArgument((A.nRows() >= A.nCols()),
                "QR decomposition by Householder Reflection applies to matrix where number of rows > number of columns");

        this.nrows = A.nRows();
        this.ncols = A.nCols();
        this.epsilon = epsilon;

        R = new UpperTriangularMatrix(ncols);
        R.set(1, 1, 0);//allocate space
        cols = new Vector[ncols + 1];//0-th element is not used
        Hs = new Householder[ncols + 1];//0-th element is not used

        /*
         * Make a copy of A by columns
         */
        for (int i = 1; i <= ncols; ++i) {
            cols[i] = A.column(i);
        }

        /*
         * We apply the Householder operator repeatedly on each columns,
         * hence looping over columns
         */
        for (int i = 1; i <= Math.min(ncols, nrows - 1); i++) {//if square, i <= nRows-1; else i <= nCols (b/c nRows > nCols by assumption)
            /*
             * We will zero out in this column all entries below the i-th entry.
             * E.g., for the first column, we want to have this form (a,0,0,..)'.
             *
             * We work on a sub-matrix., namely,
             * A[i:nRows, i:nCols]
             */
            Context defn = Householder.getContext(subVector(cols[i], i, nrows));
            Vector u = defn.generator;
            /*
             * Treat the (sub-)column as 0, if the norm is too small.
             * This is to prevent spurious Householder operation on a "zero" column,
             * which may produce incorrect result.
             */
            if (compare(Math.abs(defn.lambda), epsilon) <= 0) {//linear dependence on previous columns
                u = new DenseVector(nrows);//ZERO vector
                Hs[i] = new Householder(u);
                R.set(i, i, 0);//record R's diagonal
            } else {//to reduce computation, we perform Householder reflection only for linearly independent columns
                u = concat(
                        new DenseVector(i - 1),//ZERO vector
                        u);
                Hs[i] = new Householder(u);

                /*
                 * Apply the householder transformation to the rest of the columns, but only partially (sub-matrix).
                 * That is, we set only a sub-matrix.
                 *
                 * TODO: how to improve performance by doing inner-product (reflect) only on non-zero entries?
                 */
                for (int j = i + 1; j <= ncols; ++j) {
                    cols[j] = Hs[i].reflect(cols[j]);
                }

                R.set(i, i, defn.lambda);//record R's diagonal
            }
        }

        computeR();
    }

    /**
     * Construct an instance of the Householder reflection process to orthogonalize a matrix.
     *
     * @param A a matrix, where #rows >= #columns
     */
    public HouseholderReflection(Matrix A) {
        this(A, SuanShuUtils.autoEpsilon(DenseMatrix.as(A).to1DArray()));
    }

    /**
     * Compute <code>R</code>.
     */
    private void computeR() {
        /*
         * R's diagonal is already correct by the constructor.
         * Just copy the upper triangle entries from the Cols,
         * which are A modified by a sequence of Householder transformations.
         */
        for (int i = 2; i <= ncols; ++i) {//loop over columns
            for (int j = 1; j < i; ++j) {//loop over rows
                R.set(j, i, cols[i].get(j));
            }
        }

        //for square matrix, last column was not worked on to generate diagonal element
        if (nrows == ncols) {
            R.set(ncols, ncols, cols[ncols].get(ncols));
        }
    }

    /**
     * Get a copy of the <code>Q</code> matrix in the QR decomposition.
     *
     * <blockquote><code><pre>
     * Qn %*% ... %*% Q3 %*% Q2 %*% Q1 %*% A = R
     * A = [Qn %*% ... %*% Q3 %*% Q2 %*% Q1]^-1 %*% R
     * A = [Q1^-1 %*% Q2^-1 %*% ... %*% Qn^-1] %*% R
     * A = [Q1 %*% Q2 %*% ... %*% Qn] %*% R, by the unitary property of Householder matrix
     * So, Q = Q1 %*% Q2 %*% ... %*% Qn.
     * </pre></code></blockquote>
     * 
     * <p>
     * To compute <code>Q</code>,
     * instead of explicitly doing this multiplication, we can improve the performance
     * by applying <code>Q<sub>i</sub></code>'s repeatedly on an identity matrix.
     *
     * <p>
     * Dimension of <code>Q</code> is <code>nRows x nCols</code>.
     *
     * @return a copy of the <code>Q</code> matrix in the QR decomposition
     */
    @Override
    public Matrix Q() {
        Matrix Q = Householder.product(Hs, 1, ncols, nrows, ncols);//TODO: cache Q
        return Q;
    }

    @Override
    public UpperTriangularMatrix R() {
        return new UpperTriangularMatrix(R);
    }

    /**
     * Get a copy of the <code>P</code> pivoting matrix in the QR decomposition.
     * <p>
     * Householder process needs no pivoting.
     * Hence, <code>P</code> is always an identity matrix.
     *
     * @return an identity matrix
     */
    @Override
    public PermutationMatrix P() {
        return new PermutationMatrix(ncols);
    }

    /**
     * Rank is computed by counting the number of non-zero rows in <code>R</code>.
     *
     * <p>
     * <code>Q</code> may have more columns than the rank.
     * The user is encouraged to interpret the result with caution.
     *
     * @return the rank
     */
    @Override
    public int rank() {
        //count the number of non-zero rows in R
        //i.e., count the number of non-zero diagonal elements in R
        int result = 0;
        for (int i = 1; i <= R.nCols(); ++i) {
            if (compare(R.get(i, i), 0, epsilon) != 0) {
                ++result;
            }
        }
        return result;
    }

    @Override
    public Matrix squareQ() {
        Matrix squareQ = Householder.product(Hs, 1, ncols, nrows, nrows);//the dimension is different from that in Q
        return squareQ;
    }

    @Override
    public Matrix tallR() {
        Matrix tallR = new DenseMatrix(nrows, ncols).ZERO();
        CreateMatrix.replace(tallR, 1, ncols, 1, ncols, R);
        return tallR;
    }
}