home | career | drupal | java | mac | mysql | perl | scala | uml | unix  

Java example source code file (BiDiagonalTransformer.java)

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

Learn more about this Java project at its project page.

Java - Java tags/keywords

bidiagonaltransformer, realmatrix

The BiDiagonalTransformer.java Java example 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.math3.linear;

import org.apache.commons.math3.util.FastMath;


/**
 * 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>
 * @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.
     */
    BiDiagonalTransformer(RealMatrix matrix) {

        final int m = matrix.getRowDimension();
        final int n = matrix.getColumnDimension();
        final int p = FastMath.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; double[][] ua = new double[m][m]; // fill up the part of the matrix not affected by Householder transforms for (int k = m - 1; k >= p; --k) { ua[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]; ua[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 -= ua[i][j] * householderVectors[i][k - diagOffset]; } alpha /= diagonal[k - diagOffset] * hK[k - diagOffset]; for (int i = k; i < m; ++i) { ua[i][j] += -alpha * householderVectors[i][k - diagOffset]; } } } } if (diagOffset > 0) { ua[0][0] = 1; } cachedU = MatrixUtils.createRealMatrix(ua); } // 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; double[][] ba = new double[m][n]; for (int i = 0; i < main.length; ++i) { ba[i][i] = main[i]; if (m < n) { if (i > 0) { ba[i][i-1] = secondary[i - 1]; } } else { if (i < main.length - 1) { ba[i][i+1] = secondary[i]; } } } cachedB = MatrixUtils.createRealMatrix(ba); } // 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; double[][] va = new double[n][n]; // fill up the part of the matrix not affected by Householder transforms for (int k = n - 1; k >= p; --k) { va[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]; va[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 -= va[i][j] * hK[i]; } beta /= diagonal[k - diagOffset] * hK[k]; for (int i = k; i < n; ++i) { va[i][j] += -beta * hK[i]; } } } } if (diagOffset > 0) { va[0][0] = 1; } cachedV = MatrixUtils.createRealMatrix(va); } // 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) ? -FastMath.sqrt(xNormSqr) : FastMath.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) ? -FastMath.sqrt(xNormSqr) : FastMath.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) ? -FastMath.sqrt(xNormSqr) : FastMath.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) ? -FastMath.sqrt(xNormSqr) : FastMath.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 Java examples (source code examples)

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



my book on functional programming

 

new blog posts

 

Copyright 1998-2019 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.