diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/linear/BiDiagonalTransformer.java')
-rw-r--r-- | src/main/java/org/apache/commons/math/linear/BiDiagonalTransformer.java | 381 |
1 files changed, 381 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/linear/BiDiagonalTransformer.java b/src/main/java/org/apache/commons/math/linear/BiDiagonalTransformer.java new file mode 100644 index 0000000..78230de --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/BiDiagonalTransformer.java @@ -0,0 +1,381 @@ +/* + * 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; + +import org.apache.commons.math.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> + * <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: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + * @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 = 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.</p> + * @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.</p> + * @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) ? -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]; + } + } + } + } + + } + } + +} |