You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2012/05/05 19:53:32 UTC
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
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
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
Posted by Luc Maisonobe <Lu...@free.fr>.
Le 05/05/2012 20:08, Thomas Neidhart a écrit :
> On 05/05/2012 08:02 PM, Luc Maisonobe wrote:
>> 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.
>
> The other parts are also working:
>
> * SchurTransformer
> * Final Eigenvalue / vector computation
>
> It needs a bit more cleanup but I plan to commit it till tomorrow.
>
> Where I am still struggling is with the internal DecompositionSolver,
> which right now only works for real eigenvalues. In case of complex
> eigenvalues, this would need to change and may result in a complex result.
>
> So I do not know if we should throw a NotSupportedException in case of
> complex eigenvalues when calling EigenDecomposition.getSolver()?
I don't think so. We could return the real eigenvalues (if there are
some), or we could return something like illegal state.
Luc
>
> Thomas
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
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
Posted by Thomas Neidhart <th...@gmail.com>.
On 05/05/2012 08:02 PM, Luc Maisonobe wrote:
> 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.
The other parts are also working:
* SchurTransformer
* Final Eigenvalue / vector computation
It needs a bit more cleanup but I plan to commit it till tomorrow.
Where I am still struggling is with the internal DecompositionSolver,
which right now only works for real eigenvalues. In case of complex
eigenvalues, this would need to change and may result in a complex result.
So I do not know if we should throw a NotSupportedException in case of
complex eigenvalues when calling EigenDecomposition.getSolver()?
Thomas
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
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
Posted by Luc Maisonobe <Lu...@free.fr>.
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 × 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
>
>
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org