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