commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From t.@apache.org
Subject 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 17:53:32 GMT
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.

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



Mime
View raw message