commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <Luc.Maison...@free.fr>
Subject Re: svn commit: r1334463 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/linear/HessenbergTransformer.java test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
Date Sat, 05 May 2012 18:02:24 GMT
Le 05/05/2012 19:53, tn@apache.org a écrit :
> Author: tn
> Date: Sat May  5 17:53:32 2012
> New Revision: 1334463
> 
> URL: http://svn.apache.org/viewvc?rev=1334463&view=rev
> Log:
> [MATH-235] add HessenbergTransformer to transform a general real matrix to Hessenberg
form.
> This is a first step for the eigenvalue decomposition of a non-symmetric matrix.

Great!

Thanks, Thomas.

Luc

> 
> Added:
>     commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java
  (with props)
>     commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
  (with props)
> 
> Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java
> URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java?rev=1334463&view=auto
> ==============================================================================
> --- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java
(added)
> +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java
Sat May  5 17:53:32 2012
> @@ -0,0 +1,233 @@
> +/*
> + * 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 org.apache.commons.math3.util.FastMath;
> +import org.apache.commons.math3.util.Precision;
> +
> +/**
> + * Class transforming a general real matrix to Hessenberg form.
> + * <p>A m &times; m matrix A can be written as the product of three matrices:
A = P
> + * &times; H &times; P<sup>T</sup> with P an unitary matrix and
H a Hessenberg
> + * matrix. Both P and H are m &times; m matrices.</p>
> + * <p>Transformation to Hessenberg form 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>
> + * <p>This class is based on the method orthes in class EigenvalueDecomposition
> + * from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
> + *
> + * @see <a href="http://mathworld.wolfram.com/HessenbergDecomposition.html">MathWorld</a>
> + * @see <a href="http://en.wikipedia.org/wiki/Householder_transformation">Householder
Transformations</a>
> + * @version $Id$
> + * @since 3.1
> + */
> +class HessenbergTransformer {
> +    /** Householder vectors. */
> +    private final double householderVectors[][];
> +    /** Temporary storage vector. */
> +    private final double ort[];
> +    /** Cached value of P. */
> +    private RealMatrix cachedP;
> +    /** Cached value of Pt. */
> +    private RealMatrix cachedPt;
> +    /** Cached value of H. */
> +    private RealMatrix cachedH;
> +
> +    /**
> +     * Build the transformation to Hessenberg form of a general matrix.
> +     *
> +     * @param matrix matrix to transform.
> +     * @throws NonSquareMatrixException if the matrix is not square.
> +     */
> +    public HessenbergTransformer(RealMatrix matrix) {
> +        if (!matrix.isSquare()) {
> +            throw new NonSquareMatrixException(matrix.getRowDimension(),
> +                    matrix.getColumnDimension());
> +        }
> +
> +        final int m = matrix.getRowDimension();
> +        householderVectors = matrix.getData();
> +        ort = new double[m];
> +        cachedP = null;
> +        cachedPt = null;
> +        cachedH = null;
> +
> +        // transform matrix
> +        transform();
> +    }
> +
> +    /**
> +     * Returns the matrix P of the transform.
> +     * <p>P is an unitary matrix, i.e. its inverse is also its transpose.</p>
> +     *
> +     * @return the P matrix
> +     */
> +    public RealMatrix getP() {
> +        if (cachedP == null) {
> +            final int n = householderVectors.length;
> +            final int high = n - 1;
> +            final double[][] pa = new double[n][n];
> +
> +            for (int i = 0; i < n; i++) {
> +                for (int j = 0; j < n; j++) {
> +                    pa[i][j] = (i == j) ? 1 : 0;
> +                }
> +            }
> +
> +            for (int m = high - 1; m >= 1; m--) {
> +                if (householderVectors[m][m - 1] != 0.0) {
> +                    for (int i = m + 1; i <= high; i++) {
> +                        ort[i] = householderVectors[i][m - 1];
> +                    }
> +
> +                    for (int j = m; j <= high; j++) {
> +                        double g = 0.0;
> +
> +                        for (int i = m; i <= high; i++) {
> +                            g += ort[i] * pa[i][j];
> +                        }
> +
> +                        // Double division avoids possible underflow
> +                        g = (g / ort[m]) / householderVectors[m][m - 1];
> +
> +                        for (int i = m; i <= high; i++) {
> +                            pa[i][j] += g * ort[i];
> +                        }
> +                    }
> +                }
> +            }
> +
> +            cachedP = MatrixUtils.createRealMatrix(pa);
> +        }
> +        return cachedP;
> +    }
> +
> +    /**
> +     * Returns the transpose of the matrix P of the transform.
> +     * <p>P is an unitary matrix, i.e. its inverse is also its transpose.</p>
> +     *
> +     * @return the transpose of the P matrix
> +     */
> +    public RealMatrix getPT() {
> +        if (cachedPt == null) {
> +            cachedPt = getP().transpose();
> +        }
> +
> +        // return the cached matrix
> +        return cachedPt;
> +    }
> +
> +    /**
> +     * Returns the Hessenberg matrix H of the transform.
> +     *
> +     * @return the H matrix
> +     */
> +    public RealMatrix getH() {
> +        if (cachedH == null) {
> +            final int m = householderVectors.length;
> +            final double[][] h = new double[m][m];
> +            for (int i = 0; i < m; ++i) {
> +                if (i > 0) {
> +                    // copy the entry of the lower sub-diagonal
> +                    h[i][i - 1] = householderVectors[i][i - 1];
> +                }
> +
> +                // copy upper triangular part of the matrix
> +                for (int j = i; j < m; ++j) {
> +                    h[i][j] = householderVectors[i][j];
> +                }
> +            }
> +            cachedH = MatrixUtils.createRealMatrix(h);
> +        }
> +
> +        // return the cached matrix
> +        return cachedH;
> +    }
> +
> +    /**
> +     * 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;
> +    }
> +
> +    /**
> +     * Transform original matrix to Hessenberg form.
> +     * <p>Transformation is done using Householder transforms.</p>
> +     */
> +    private void transform() {
> +        final int n = householderVectors.length;
> +        final int high = n - 1;
> +
> +        for (int m = 1; m <= high - 1; m++) {
> +            // Scale column.
> +            double scale = 0;
> +            for (int i = m; i <= high; i++) {
> +                scale += FastMath.abs(householderVectors[i][m - 1]);
> +            }
> +
> +            if (!Precision.equals(scale, 0)) {
> +                // Compute Householder transformation.
> +                double h = 0;
> +                for (int i = high; i >= m; i--) {
> +                    ort[i] = householderVectors[i][m - 1] / scale;
> +                    h += ort[i] * ort[i];
> +                }
> +                final double g = (ort[m] > 0) ? -FastMath.sqrt(h) : FastMath.sqrt(h);
> +
> +                h = h - ort[m] * g;
> +                ort[m] = ort[m] - g;
> +
> +                // Apply Householder similarity transformation
> +                // H = (I - u*u' / h) * H * (I - u*u' / h)
> +
> +                for (int j = m; j < n; j++) {
> +                    double f = 0;
> +                    for (int i = high; i >= m; i--) {
> +                        f += ort[i] * householderVectors[i][j];
> +                    }
> +                    f = f / h;
> +                    for (int i = m; i <= high; i++) {
> +                        householderVectors[i][j] -= f * ort[i];
> +                    }
> +                }
> +
> +                for (int i = 0; i <= high; i++) {
> +                    double f = 0;
> +                    for (int j = high; j >= m; j--) {
> +                        f += ort[j] * householderVectors[i][j];
> +                    }
> +                    f = f / h;
> +                    for (int j = m; j <= high; j++) {
> +                        householderVectors[i][j] -= f * ort[j];
> +                    }
> +                }
> +
> +                ort[m] = scale * ort[m];
> +                householderVectors[m][m - 1] = scale * g;
> +            }
> +        }
> +    }
> +}
> 
> Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java
> ------------------------------------------------------------------------------
>     svn:eol-style = native
> 
> Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java
> ------------------------------------------------------------------------------
>     svn:keywords = Id Revision HeadURL
> 
> Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/HessenbergTransformer.java
> ------------------------------------------------------------------------------
>     svn:mime-type = text/plain
> 
> Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
> URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java?rev=1334463&view=auto
> ==============================================================================
> --- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
(added)
> +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
Sat May  5 17:53:32 2012
> @@ -0,0 +1,160 @@
> +/*
> + * 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 org.junit.Test;
> +import org.junit.Assert;
> +
> +public class HessenbergTransformerTest {
> +
> +    private double[][] testSquare5 = {
> +            { 5, 4, 3, 2, 1 },
> +            { 1, 4, 0, 3, 3 },
> +            { 2, 0, 3, 0, 0 },
> +            { 3, 2, 1, 2, 5 },
> +            { 4, 2, 1, 4, 1 }
> +    };
> +
> +    private double[][] testSquare3 = {
> +            {  2, -1, 1 },
> +            { -1,  2, 1 },
> +            {  1, -1, 2 }
> +    };
> +
> +    // from http://eigen.tuxfamily.org/dox/classEigen_1_1HessenbergDecomposition.html
> +
> +    private double[][] testRandom = {
> +            {  0.680,  0.823, -0.4440, -0.2700 },
> +            { -0.211, -0.605,  0.1080,  0.0268 },
> +            {  0.566, -0.330, -0.0452,  0.9040 },
> +            {  0.597,  0.536,  0.2580,  0.8320 }
> +    };
> +
> +    @Test
> +    public void testNonSquare() {
> +        try {
> +            new HessenbergTransformer(MatrixUtils.createRealMatrix(new double[3][2]));
> +            Assert.fail("an exception should have been thrown");
> +        } catch (NonSquareMatrixException ime) {
> +            // expected behavior
> +        }
> +    }
> +
> +    @Test
> +    public void testAEqualPHPt() {
> +        checkAEqualPHPt(MatrixUtils.createRealMatrix(testSquare5));
> +        checkAEqualPHPt(MatrixUtils.createRealMatrix(testSquare3));
> +        checkAEqualPHPt(MatrixUtils.createRealMatrix(testRandom));
> +   }
> +
> +    private void checkAEqualPHPt(RealMatrix matrix) {
> +        HessenbergTransformer transformer = new HessenbergTransformer(matrix);
> +        RealMatrix p  = transformer.getP();
> +        RealMatrix pT = transformer.getPT();
> +        RealMatrix h  = transformer.getH();
> +        double norm = p.multiply(h).multiply(pT).subtract(matrix).getNorm();
> +        Assert.assertEquals(0, norm, 4.0e-14);
> +    }
> +
> +    @Test
> +    public void testPOrthogonal() {
> +        checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getP());
> +        checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getP());
> +    }
> +
> +    @Test
> +    public void testPTOrthogonal() {
> +        checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getPT());
> +        checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getPT());
> +    }
> +
> +    private void checkOrthogonal(RealMatrix m) {
> +        RealMatrix mTm = m.transpose().multiply(m);
> +        RealMatrix id  = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
> +        Assert.assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);
> +    }
> +
> +    @Test
> +    public void testHessenbergForm() {
> +        checkHessenbergForm(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getH());
> +        checkHessenbergForm(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getH());
> +    }
> +
> +    private void checkHessenbergForm(RealMatrix m) {
> +        final int rows = m.getRowDimension();
> +        final int cols = m.getColumnDimension();
> +        for (int i = 0; i < rows; ++i) {
> +            for (int j = 0; j < cols; ++j) {
> +                if (i > j + 1) {
> +                    Assert.assertEquals(0, m.getEntry(i, j), 1.0e-16);
> +                }
> +            }
> +        }
> +    }
> +
> +    @Test
> +    public void testMatricesValues5() {
> +        checkMatricesValues(testSquare5,
> +                            new double[][] {
> +                                { 1.0,  0.0,                0.0,                0.0,
               0.0               },
> +                                { 0.0, -0.182574185835055,  0.784218758628863,  0.395029040913988,
-0.442289115981669 },
> +                                { 0.0, -0.365148371670111, -0.337950625265477, -0.374110794088820,
-0.782621974707823 },
> +                                { 0.0, -0.547722557505166,  0.402941130124223, -0.626468266309003,
 0.381019628053472 },
> +                                { 0.0, -0.730296743340221, -0.329285224617644,  0.558149336547665,
 0.216118545309225 }
> +                            },
> +                            new double[][] {
> +                                {  5.0,              -3.65148371670111,  2.59962019434982,
-0.237003414680848, -3.13886458663398  },
> +                                { -5.47722557505166,  6.9,              -2.29164066120599,
 0.207283564429169,  0.703858369151728 },
> +                                {  0.0,              -4.21386600008432,  2.30555659846067,
 2.74935928725112,   0.857569835914113 },
> +                                {  0.0,               0.0,               2.86406180891882,
-1.11582249161595,   0.817995267184158 },
> +                                {  0.0,               0.0,               0.0,      
        0.683518597386085,  1.91026589315528  }
> +                            });
> +    }
> +
> +    @Test
> +    public void testMatricesValues3() {
> +        checkMatricesValues(testSquare3,
> +                            new double[][] {
> +                                {  1.0,  0.0,               0.0               },
> +                                {  0.0, -0.707106781186547, 0.707106781186547 },
> +                                {  0.0,  0.707106781186547, 0.707106781186548 },
> +                            },
> +                            new double[][] {
> +                                {  2.0,              1.41421356237309,  0.0 },
> +                                {  1.41421356237310, 2.0,              -1.0 },
> +                                {  0.0,              1.0,               2.0 },
> +                            });
> +    }
> +
> +    private void checkMatricesValues(double[][] matrix, double[][] pRef, double[][]
hRef) {
> +
> +        HessenbergTransformer transformer =
> +            new HessenbergTransformer(MatrixUtils.createRealMatrix(matrix));
> +
> +        // check values against known references
> +        RealMatrix p = transformer.getP();
> +        Assert.assertEquals(0, p.subtract(MatrixUtils.createRealMatrix(pRef)).getNorm(),
1.0e-14);
> +
> +        RealMatrix h = transformer.getH();
> +        Assert.assertEquals(0, h.subtract(MatrixUtils.createRealMatrix(hRef)).getNorm(),
1.0e-14);
> +
> +        // check the same cached instance is returned the second time
> +        Assert.assertTrue(p == transformer.getP());
> +        Assert.assertTrue(h == transformer.getH());
> +    }
> +}
> 
> Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
> ------------------------------------------------------------------------------
>     svn:eol-style = native
> 
> Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
> ------------------------------------------------------------------------------
>     svn:keywords = Id Revision HeadURL
> 
> Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/HessenbergTransformerTest.java
> ------------------------------------------------------------------------------
>     svn:mime-type = text/plain
> 
> 
> 


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message