All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.apache.commons.math3.linear.BiDiagonalTransformer Maven / Gradle / Ivy

Go to download

The Apache Commons Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang.

There is a newer version: 62
Show newest version
/*
 * 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.
 * 

Any m × n matrix A can be written as the product of three matrices: * A = U × B × VT 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.

*

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.

* @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. *

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. *

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. *

Note that since this class is only intended for internal use, * it returns directly a reference to its internal arrays, not a copy.

* @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. *

Note that since this class is only intended for internal use, * it returns directly a reference to its internal arrays, not a copy.

* @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. *

Note that since this class is only intended for internal use, * it returns directly a reference to its internal arrays, not a copy.

* @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. *

Transformation is done using alternate Householder transforms * on columns and rows.

*/ 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. *

Transformation is done using alternate Householder transforms * on rows and columns.

*/ 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]; } } } } } } }




© 2015 - 2024 Weber Informatics LLC | Privacy Policy