diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java | 270 |
1 files changed, 270 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java b/src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java new file mode 100644 index 0000000..9260aa0 --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java @@ -0,0 +1,270 @@ +/* + * 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 java.util.Arrays; + +import org.apache.commons.math.util.FastMath; + + +/** + * Class transforming a symmetrical matrix to tridiagonal shape. + * <p>A symmetrical m × m matrix A can be written as the product of three matrices: + * A = Q × T × Q<sup>T</sup> with Q an orthogonal matrix and T a symmetrical + * tridiagonal matrix. Both Q and T are m × m matrices.</p> + * <p>This implementation only uses the upper part of the matrix, the part below the + * diagonal is not accessed at all.</p> + * <p>Transformation to tridiagonal shape is often not a goal by itself, but it is + * an intermediate step in more general decomposition algorithms like {@link + * EigenDecomposition eigen 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 TriDiagonalTransformer { + + /** Householder vectors. */ + private final double householderVectors[][]; + + /** Main diagonal. */ + private final double[] main; + + /** Secondary diagonal. */ + private final double[] secondary; + + /** Cached value of Q. */ + private RealMatrix cachedQ; + + /** Cached value of Qt. */ + private RealMatrix cachedQt; + + /** Cached value of T. */ + private RealMatrix cachedT; + + /** + * Build the transformation to tridiagonal shape of a symmetrical matrix. + * <p>The specified matrix is assumed to be symmetrical without any check. + * Only the upper triangular part of the matrix is used.</p> + * @param matrix the symmetrical matrix to transform. + * @exception InvalidMatrixException if matrix is not square + */ + public TriDiagonalTransformer(RealMatrix matrix) + throws InvalidMatrixException { + if (!matrix.isSquare()) { + throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension()); + } + + final int m = matrix.getRowDimension(); + householderVectors = matrix.getData(); + main = new double[m]; + secondary = new double[m - 1]; + cachedQ = null; + cachedQt = null; + cachedT = null; + + // transform matrix + transform(); + + } + + /** + * Returns the matrix Q of the transform. + * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p> + * @return the Q matrix + */ + public RealMatrix getQ() { + if (cachedQ == null) { + cachedQ = getQT().transpose(); + } + return cachedQ; + } + + /** + * Returns the transpose of the matrix Q of the transform. + * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p> + * @return the Q matrix + */ + public RealMatrix getQT() { + + if (cachedQt == null) { + + final int m = householderVectors.length; + cachedQt = MatrixUtils.createRealMatrix(m, m); + + // build up first part of the matrix by applying Householder transforms + for (int k = m - 1; k >= 1; --k) { + final double[] hK = householderVectors[k - 1]; + final double inv = 1.0 / (secondary[k - 1] * hK[k]); + cachedQt.setEntry(k, k, 1); + if (hK[k] != 0.0) { + double beta = 1.0 / secondary[k - 1]; + cachedQt.setEntry(k, k, 1 + beta * hK[k]); + for (int i = k + 1; i < m; ++i) { + cachedQt.setEntry(k, i, beta * hK[i]); + } + for (int j = k + 1; j < m; ++j) { + beta = 0; + for (int i = k + 1; i < m; ++i) { + beta += cachedQt.getEntry(j, i) * hK[i]; + } + beta *= inv; + cachedQt.setEntry(j, k, beta * hK[k]); + for (int i = k + 1; i < m; ++i) { + cachedQt.addToEntry(j, i, beta * hK[i]); + } + } + } + } + cachedQt.setEntry(0, 0, 1); + + } + + // return the cached matrix + return cachedQt; + + } + + /** + * Returns the tridiagonal matrix T of the transform. + * @return the T matrix + */ + public RealMatrix getT() { + + if (cachedT == null) { + + final int m = main.length; + cachedT = MatrixUtils.createRealMatrix(m, m); + for (int i = 0; i < m; ++i) { + cachedT.setEntry(i, i, main[i]); + if (i > 0) { + cachedT.setEntry(i, i - 1, secondary[i - 1]); + } + if (i < main.length - 1) { + cachedT.setEntry(i, i + 1, secondary[i]); + } + } + + } + + // return the cached matrix + return cachedT; + + } + + /** + * 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 T 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 T matrix + */ + double[] getMainDiagonalRef() { + return main; + } + + /** + * Get the secondary diagonal elements of the matrix T 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 T matrix + */ + double[] getSecondaryDiagonalRef() { + return secondary; + } + + /** + * Transform original matrix to tridiagonal form. + * <p>Transformation is done using Householder transforms.</p> + */ + private void transform() { + + final int m = householderVectors.length; + final double[] z = new double[m]; + for (int k = 0; k < m - 1; k++) { + + //zero-out a row and a column simultaneously + final double[] hK = householderVectors[k]; + main[k] = hK[k]; + double xNormSqr = 0; + for (int j = k + 1; j < m; ++j) { + final double c = hK[j]; + xNormSqr += c * c; + } + final double a = (hK[k + 1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); + secondary[k] = a; + if (a != 0.0) { + // apply Householder transform from left and right simultaneously + + hK[k + 1] -= a; + final double beta = -1 / (a * hK[k + 1]); + + // compute a = beta A v, where v is the Householder vector + // this loop is written in such a way + // 1) only the upper triangular part of the matrix is accessed + // 2) access is cache-friendly for a matrix stored in rows + Arrays.fill(z, k + 1, m, 0); + for (int i = k + 1; i < m; ++i) { + final double[] hI = householderVectors[i]; + final double hKI = hK[i]; + double zI = hI[i] * hKI; + for (int j = i + 1; j < m; ++j) { + final double hIJ = hI[j]; + zI += hIJ * hK[j]; + z[j] += hIJ * hKI; + } + z[i] = beta * (z[i] + zI); + } + + // compute gamma = beta vT z / 2 + double gamma = 0; + for (int i = k + 1; i < m; ++i) { + gamma += z[i] * hK[i]; + } + gamma *= beta / 2; + + // compute z = z - gamma v + for (int i = k + 1; i < m; ++i) { + z[i] -= gamma * hK[i]; + } + + // update matrix: A = A - v zT - z vT + // only the upper triangular part of the matrix is updated + for (int i = k + 1; i < m; ++i) { + final double[] hI = householderVectors[i]; + for (int j = i; j < m; ++j) { + hI[j] -= hK[i] * z[j] + z[i] * hK[j]; + } + } + + } + + } + main[m - 1] = householderVectors[m - 1][m - 1]; + } + +} |