diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java | 388 |
1 files changed, 388 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java b/src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java new file mode 100644 index 0000000..063b202 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java @@ -0,0 +1,388 @@ +/* + * 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</sup> 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>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. + * + * <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. + * + * @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. + * + * @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. + * + * @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. + */ + 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. + */ + 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]; + } + } + } + } + } + } +} |