commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/TriDiagonalTransformer.java - commons-math - Git at Google

 /*
  * 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.math4.legacy.linear;

 import java.util.Arrays;

 import org.apache.commons.math4.core.jdkmath.JdkMath;


 /**
  * Class transforming a symmetrical matrix to tridiagonal shape.
  * <p>A symmetrical m &times; m matrix A can be written as the product of three matrices:
  * A = Q &times; T &times; Q<sup>T</sup> with Q an orthogonal matrix and T a symmetrical
  * tridiagonal matrix. Both Q and T are m &times; 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.</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;
             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.</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) ? -JdkMath.sqrt(xNormSqr) : JdkMath.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];
     }
 }
	/*
	* 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.math4.legacy.linear;

	import java.util.Arrays;

	import org.apache.commons.math4.core.jdkmath.JdkMath;


	/**
	* 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>
	* @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.</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;
	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.</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) ? -JdkMath.sqrt(xNormSqr) : JdkMath.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];
	}
	}