home | career | drupal | java | mac | mysql | perl | scala | uml | unix  

Java example source code file (TriDiagonalTransformer.java)

This example Java source code file (TriDiagonalTransformer.java) is included in the alvinalexander.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

Learn more about this Java project at its project page.

Java - Java tags/keywords

realmatrix, tridiagonaltransformer, util

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/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.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++) { //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]; } }

Other Java examples (source code examples)

Here is a short list of links related to this Java TriDiagonalTransformer.java source code file:



my book on functional programming

 

new blog posts

 

Copyright 1998-2019 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.