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