org.ejml.alg.dense.decomposition.bidiagonal.BidiagonalDecompositionNaive Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of ejml Show documentation
A fast and easy to use dense matrix linear algebra library written in Java.
There is a newer version: 0.25
/*
 * Copyright (c) 2009-2012, Peter Abeles. All Rights Reserved.
 *
 * This file is part of Efficient Java Matrix Library (EJML).
 *
 * EJML is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as
 * published by the Free Software Foundation, either version 3
 * of the License, or (at your option) any later version.
 *
 * EJML is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with EJML.  If not, see .
 */

package org.ejml.alg.dense.decomposition.bidiagonal;

import org.ejml.data.DenseMatrix64F;
import org.ejml.ops.SpecializedOps;
import org.ejml.simple.SimpleMatrix;


/**
 * A slower but much simpler version of {@link BidiagonalDecompositionRow} that internally uses
 * SimpleMatrix and explicitly computes the householder matrices.  This was easier to code up and is
 * used to validate other implementations.
 *
 * @author Peter Abeles
 */
public class BidiagonalDecompositionNaive {
    private SimpleMatrix U;
    private SimpleMatrix B;
    private SimpleMatrix V;

    // number of rows
    private int m;
    // number of columns
    private int n;
    // smallest of m and n
    private int min;


    DenseMatrix64F u;

    public SimpleMatrix getU() {
        return U;
    }

    public SimpleMatrix getB() {
        return B;
    }

    public SimpleMatrix getV() {
        return V;
    }

    /**
     * Computes the decomposition of the provided matrix.  If no errors are detected then true is returned,
     * false otherwise.
     * @param A  The matrix that is being decomposed.  Not modified.
     * @return If it detects any errors or not.
     */
    public boolean decompose( DenseMatrix64F A )
    {
        init(A);
        return _decompose();
    }

    protected void init(DenseMatrix64F A) {
        m = A.numRows;
        n = A.numCols;

        min = Math.min(m,n);

        U = SimpleMatrix.identity(m);
        B = new SimpleMatrix(A);
        V = SimpleMatrix.identity(n);

        int max = Math.max(m,n);
        u = new DenseMatrix64F(max,1);
    }



    /**
     * Internal function for computing the decomposition.
     */
    private boolean _decompose() {
        for( int k = 0; k < min; k++ ) {
            computeU(k);
            computeV(k);
        }

        return true;
    }

    protected void computeU( int k) {
        u.reshape(m,1, false);
        double u[] = this.u.data;

        // find the largest value in this column
        // this is used to normalize the column and mitigate overflow/underflow
        double max = 0;

        for( int i = k; i < m; i++ ) {
            // copy the householder vector to vector outside of the matrix to reduce caching issues
            // big improvement on larger matrices and a relatively small performance hit on small matrices.
            double val = u[i] = B.get(i,k);
            val = Math.abs(val);
            if( val > max )
                max = val;
        }

        if( max > 0 ) {
            // -------- set up the reflector Q_k

            double tau = 0;
            // normalize to reduce overflow/underflow
            // and compute tau for the reflector
            for( int i = k; i < m; i++ ) {
                double val = u[i] /= max;
                tau += val*val;
            }

            tau = Math.sqrt(tau);

            if( u[k] < 0 )
                tau = -tau;

            // write the reflector into the lower left column of the matrix
            double nu = u[k] + tau;
            u[k] = 1.0;

            for( int i = k+1; i < m; i++ ) {
                u[i] /= nu;
            }

            SimpleMatrix Q_k = SimpleMatrix.wrap(SpecializedOps.createReflector(this.u,nu/tau));
            U = U.mult(Q_k);
            B = Q_k.mult(B);
        }
    }

    protected void computeV(int k) {
        u.reshape(n,1, false);
        u.zero();
        double u[] = this.u.data;


        // find the largest value in this column
        // this is used to normalize the column and mitigate overflow/underflow
        double max = 0;
        
        for( int i = k+1; i < n; i++ ) {
            // copy the householder vector to vector outside of the matrix to reduce caching issues
            // big improvement on larger matrices and a relatively small performance hit on small matrices.
            double val = u[i] = B.get(k,i);
            val = Math.abs(val);
            if( val > max )
                max = val;
        }

        if( max > 0 ) {
            // -------- set up the reflector Q_k

            double tau = 0;
            // normalize to reduce overflow/underflow
            // and compute tau for the reflector
            for( int i = k+1; i < n; i++ ) {
                double val = u[i] /= max;
                tau += val*val;
            }

            tau = Math.sqrt(tau);

            if( u[k+1] < 0 )
                tau = -tau;

            // write the reflector into the lower left column of the matrix
            double nu = u[k+1] + tau;
            u[k+1] = 1.0;

            for( int i = k+2; i < n; i++ ) {
                u[i] /= nu;
            }

            // ---------- multiply on the left by Q_k
            SimpleMatrix Q_k = SimpleMatrix.wrap(SpecializedOps.createReflector(this.u,nu/tau));

            V = V.mult(Q_k);
            B = B.mult(Q_k);
        }
    }

}