

Java example source code file (TriDiagonalTransformer.java)
The TriDiagonalTransformer.java Java example source code/* * 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/LICENSE2.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 java.util.Arrays; import org.apache.commons.math3.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 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> * @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 Symmetrical matrix to transform. * @throws NonSquareMatrixException if the matrix is not square. */ TriDiagonalTransformer(RealMatrix matrix) { 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. * @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. * @return the Q matrix */ public RealMatrix getQT() { if (cachedQt == null) { final int m = householderVectors.length; double[][] qta = new double[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]; qta[k][k] = 1; if (hK[k] != 0.0) { final double inv = 1.0 / (secondary[k  1] * hK[k]); double beta = 1.0 / secondary[k  1]; qta[k][k] = 1 + beta * hK[k]; for (int i = k + 1; i < m; ++i) { qta[k][i] = beta * hK[i]; } for (int j = k + 1; j < m; ++j) { beta = 0; for (int i = k + 1; i < m; ++i) { beta += qta[j][i] * hK[i]; } beta *= inv; qta[j][k] = beta * hK[k]; for (int i = k + 1; i < m; ++i) { qta[j][i] += beta * hK[i]; } } } } qta[0][0] = 1; cachedQt = MatrixUtils.createRealMatrix(qta); } // 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; double[][] ta = new double[m][m]; for (int i = 0; i < m; ++i) { ta[i][i] = main[i]; if (i > 0) { ta[i][i  1] = secondary[i  1]; } if (i < main.length  1) { ta[i][i + 1] = secondary[i]; } } cachedT = MatrixUtils.createRealMatrix(ta); } // 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. */ private void transform() { final int m = householderVectors.length; final double[] z = new double[m]; for (int k = 0; k < m  1; k++) { //zeroout 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 cachefriendly 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]; } } Other Java examples (source code examples)Here is a short list of links related to this Java TriDiagonalTransformer.java source code file: 
Copyright 19982019 Alvin Alexander, alvinalexander.com
All Rights Reserved.
A percentage of advertising revenue from
pages under the /java/jwarehouse
URI on this website is
paid back to open source projects.