Le 05/05/2012 19:53, t...@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-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
