alvinalexander.com | career | drupal | java | mac | mysql | perl | scala | uml | unix  

Commons Math example source code file (BiDiagonalTransformer.java)

This example Commons Math source code file (BiDiagonalTransformer.java) is included in the DevDaily.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

Java - Commons Math tags/keywords

bidiagonaltransformer, bidiagonaltransformer, realmatrix, realmatrix

The Commons Math BiDiagonalTransformer.java source code

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You 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 org.apache.commons.math.linear;


/**
 * Class transforming any matrix to bi-diagonal shape.
 * <p>Any m × n matrix A can be written as the product of three matrices:
 * A = U × B × V<sup>T with U an m × m orthogonal matrix,
 * B an m × n bi-diagonal matrix (lower diagonal if m < n, upper diagonal
 * otherwise), and V an n × n orthogonal matrix.</p>
 * <p>Transformation to bi-diagonal shape is often not a goal by itself, but it is
 * an intermediate step in more general decomposition algorithms like {@link
 * SingularValueDecomposition Singular Value Decomposition}. This class is therefore
 * intended for internal use by the library and is not public. As a consequence of
 * this explicitly limited scope, many methods directly returns references to
 * internal arrays, not copies.</p>
 * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
 * @since 2.0
 */
class BiDiagonalTransformer {

    /** Householder vectors. */
    private final double householderVectors[][];

    /** Main diagonal. */
    private final double[] main;

    /** Secondary diagonal. */
    private final double[] secondary;

    /** Cached value of U. */
    private RealMatrix cachedU;

    /** Cached value of B. */
    private RealMatrix cachedB;

    /** Cached value of V. */
    private RealMatrix cachedV;

    /**
     * Build the transformation to bi-diagonal shape of a matrix.
     * @param matrix the matrix to transform.
     */
    public BiDiagonalTransformer(RealMatrix matrix) {

        final int m = matrix.getRowDimension();
        final int n = matrix.getColumnDimension();
        final int p = Math.min(m, n);
        householderVectors = matrix.getData();
        main      = new double[p];
        secondary = new double[p - 1];
        cachedU   = null;
        cachedB   = null;
        cachedV   = null;

        // transform matrix
        if (m >= n) {
            transformToUpperBiDiagonal();
        } else {
            transformToLowerBiDiagonal();
        }

    }

    /**
     * Returns the matrix U of the transform.
     * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.

* @return the U matrix */ public RealMatrix getU() { if (cachedU == null) { final int m = householderVectors.length; final int n = householderVectors[0].length; final int p = main.length; final int diagOffset = (m >= n) ? 0 : 1; final double[] diagonal = (m >= n) ? main : secondary; cachedU = MatrixUtils.createRealMatrix(m, m); // fill up the part of the matrix not affected by Householder transforms for (int k = m - 1; k >= p; --k) { cachedU.setEntry(k, k, 1); } // build up first part of the matrix by applying Householder transforms for (int k = p - 1; k >= diagOffset; --k) { final double[] hK = householderVectors[k]; cachedU.setEntry(k, k, 1); if (hK[k - diagOffset] != 0.0) { for (int j = k; j < m; ++j) { double alpha = 0; for (int i = k; i < m; ++i) { alpha -= cachedU.getEntry(i, j) * householderVectors[i][k - diagOffset]; } alpha /= diagonal[k - diagOffset] * hK[k - diagOffset]; for (int i = k; i < m; ++i) { cachedU.addToEntry(i, j, -alpha * householderVectors[i][k - diagOffset]); } } } } if (diagOffset > 0) { cachedU.setEntry(0, 0, 1); } } // return the cached matrix return cachedU; } /** * Returns the bi-diagonal matrix B of the transform. * @return the B matrix */ public RealMatrix getB() { if (cachedB == null) { final int m = householderVectors.length; final int n = householderVectors[0].length; cachedB = MatrixUtils.createRealMatrix(m, n); for (int i = 0; i < main.length; ++i) { cachedB.setEntry(i, i, main[i]); if (m < n) { if (i > 0) { cachedB.setEntry(i, i - 1, secondary[i - 1]); } } else { if (i < main.length - 1) { cachedB.setEntry(i, i + 1, secondary[i]); } } } } // return the cached matrix return cachedB; } /** * Returns the matrix V of the transform. * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.

* @return the V matrix */ public RealMatrix getV() { if (cachedV == null) { final int m = householderVectors.length; final int n = householderVectors[0].length; final int p = main.length; final int diagOffset = (m >= n) ? 1 : 0; final double[] diagonal = (m >= n) ? secondary : main; cachedV = MatrixUtils.createRealMatrix(n, n); // fill up the part of the matrix not affected by Householder transforms for (int k = n - 1; k >= p; --k) { cachedV.setEntry(k, k, 1); } // build up first part of the matrix by applying Householder transforms for (int k = p - 1; k >= diagOffset; --k) { final double[] hK = householderVectors[k - diagOffset]; cachedV.setEntry(k, k, 1); if (hK[k] != 0.0) { for (int j = k; j < n; ++j) { double beta = 0; for (int i = k; i < n; ++i) { beta -= cachedV.getEntry(i, j) * hK[i]; } beta /= diagonal[k - diagOffset] * hK[k]; for (int i = k; i < n; ++i) { cachedV.addToEntry(i, j, -beta * hK[i]); } } } } if (diagOffset > 0) { cachedV.setEntry(0, 0, 1); } } // return the cached matrix return cachedV; } /** * Get the Householder vectors of the transform. * <p>Note that since this class is only intended for internal use, * it returns directly a reference to its internal arrays, not a copy.</p> * @return the main diagonal elements of the B matrix */ double[][] getHouseholderVectorsRef() { return householderVectors; } /** * Get the main diagonal elements of the matrix B of the transform. * <p>Note that since this class is only intended for internal use, * it returns directly a reference to its internal arrays, not a copy.</p> * @return the main diagonal elements of the B matrix */ double[] getMainDiagonalRef() { return main; } /** * Get the secondary diagonal elements of the matrix B of the transform. * <p>Note that since this class is only intended for internal use, * it returns directly a reference to its internal arrays, not a copy.</p> * @return the secondary diagonal elements of the B matrix */ double[] getSecondaryDiagonalRef() { return secondary; } /** * Check if the matrix is transformed to upper bi-diagonal. * @return true if the matrix is transformed to upper bi-diagonal */ boolean isUpperBiDiagonal() { return householderVectors.length >= householderVectors[0].length; } /** * Transform original matrix to upper bi-diagonal form. * <p>Transformation is done using alternate Householder transforms * on columns and rows.</p> */ private void transformToUpperBiDiagonal() { final int m = householderVectors.length; final int n = householderVectors[0].length; for (int k = 0; k < n; k++) { //zero-out a column double xNormSqr = 0; for (int i = k; i < m; ++i) { final double c = householderVectors[i][k]; xNormSqr += c * c; } final double[] hK = householderVectors[k]; final double a = (hK[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr); main[k] = a; if (a != 0.0) { hK[k] -= a; for (int j = k + 1; j < n; ++j) { double alpha = 0; for (int i = k; i < m; ++i) { final double[] hI = householderVectors[i]; alpha -= hI[j] * hI[k]; } alpha /= a * householderVectors[k][k]; for (int i = k; i < m; ++i) { final double[] hI = householderVectors[i]; hI[j] -= alpha * hI[k]; } } } if (k < n - 1) { //zero-out a row xNormSqr = 0; for (int j = k + 1; j < n; ++j) { final double c = hK[j]; xNormSqr += c * c; } final double b = (hK[k + 1] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr); secondary[k] = b; if (b != 0.0) { hK[k + 1] -= b; for (int i = k + 1; i < m; ++i) { final double[] hI = householderVectors[i]; double beta = 0; for (int j = k + 1; j < n; ++j) { beta -= hI[j] * hK[j]; } beta /= b * hK[k + 1]; for (int j = k + 1; j < n; ++j) { hI[j] -= beta * hK[j]; } } } } } } /** * Transform original matrix to lower bi-diagonal form. * <p>Transformation is done using alternate Householder transforms * on rows and columns.</p> */ private void transformToLowerBiDiagonal() { final int m = householderVectors.length; final int n = householderVectors[0].length; for (int k = 0; k < m; k++) { //zero-out a row final double[] hK = householderVectors[k]; double xNormSqr = 0; for (int j = k; j < n; ++j) { final double c = hK[j]; xNormSqr += c * c; } final double a = (hK[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr); main[k] = a; if (a != 0.0) { hK[k] -= a; for (int i = k + 1; i < m; ++i) { final double[] hI = householderVectors[i]; double alpha = 0; for (int j = k; j < n; ++j) { alpha -= hI[j] * hK[j]; } alpha /= a * householderVectors[k][k]; for (int j = k; j < n; ++j) { hI[j] -= alpha * hK[j]; } } } if (k < m - 1) { //zero-out a column final double[] hKp1 = householderVectors[k + 1]; xNormSqr = 0; for (int i = k + 1; i < m; ++i) { final double c = householderVectors[i][k]; xNormSqr += c * c; } final double b = (hKp1[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr); secondary[k] = b; if (b != 0.0) { hKp1[k] -= b; for (int j = k + 1; j < n; ++j) { double beta = 0; for (int i = k + 1; i < m; ++i) { final double[] hI = householderVectors[i]; beta -= hI[j] * hI[k]; } beta /= b * hKp1[k]; for (int i = k + 1; i < m; ++i) { final double[] hI = householderVectors[i]; hI[j] -= beta * hI[k]; } } } } } } }

Other Commons Math examples (source code examples)

Here is a short list of links related to this Commons Math BiDiagonalTransformer.java source code file:

... this post is sponsored by my books ...

#1 New Release!

FP Best Seller

 

new blog posts

 

Copyright 1998-2021 Alvin Alexander, alvinalexander.com
All Rights Reserved.

A percentage of advertising revenue from
pages under the /java/jwarehouse URI on this website is
paid back to open source projects.