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 × m matrix A can be written as the product of three matrices:
A = P
+ * × H × P<sup>T</sup> with P an unitary matrix and H a Hessenberg
+ * matrix. Both P and H are m × 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
