View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.linear;
24  
25  import java.util.Random;
26  
27  import org.hipparchus.exception.LocalizedCoreFormats;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.random.RandomDataGenerator;
30  import org.junit.Assert;
31  import org.junit.Test;
32  
33  public class HessenbergTransformerTest {
34  
35      private double[][] testSquare5 = {
36              { 5, 4, 3, 2, 1 },
37              { 1, 4, 0, 3, 3 },
38              { 2, 0, 3, 0, 0 },
39              { 3, 2, 1, 2, 5 },
40              { 4, 2, 1, 4, 1 }
41      };
42  
43      private double[][] testSquare3 = {
44              {  2, -1, 1 },
45              { -1,  2, 1 },
46              {  1, -1, 2 }
47      };
48  
49      // from http://eigen.tuxfamily.org/dox/classEigen_1_1HessenbergDecomposition.html
50  
51      private double[][] testRandom = {
52              {  0.680,  0.823, -0.4440, -0.2700 },
53              { -0.211, -0.605,  0.1080,  0.0268 },
54              {  0.566, -0.330, -0.0452,  0.9040 },
55              {  0.597,  0.536,  0.2580,  0.8320 }
56      };
57  
58      @Test
59      public void testNonSquare() {
60          try {
61              new HessenbergTransformer(MatrixUtils.createRealMatrix(new double[3][2]));
62              Assert.fail("an exception should have been thrown");
63          } catch (MathIllegalArgumentException ime) {
64              Assert.assertEquals(LocalizedCoreFormats.NON_SQUARE_MATRIX, ime.getSpecifier());
65          }
66      }
67  
68      @Test
69      public void testAEqualPHPt() {
70          checkAEqualPHPt(MatrixUtils.createRealMatrix(testSquare5));
71          checkAEqualPHPt(MatrixUtils.createRealMatrix(testSquare3));
72          checkAEqualPHPt(MatrixUtils.createRealMatrix(testRandom));
73     }
74  
75      @Test
76      public void testPOrthogonal() {
77          checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getP());
78          checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getP());
79      }
80  
81      @Test
82      public void testPTOrthogonal() {
83          checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getPT());
84          checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getPT());
85      }
86  
87      @Test
88      public void testHessenbergForm() {
89          checkHessenbergForm(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getH());
90          checkHessenbergForm(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getH());
91      }
92  
93      @Test
94      public void testRandomData() {
95          for (int run = 0; run < 100; run++) {
96              Random r = new Random(System.currentTimeMillis());
97  
98              // matrix size
99              int size = r.nextInt(20) + 4;
100 
101             double[][] data = new double[size][size];
102             for (int i = 0; i < size; i++) {
103                 for (int j = 0; j < size; j++) {
104                     data[i][j] = r.nextInt(100);
105                 }
106             }
107 
108             RealMatrix m = MatrixUtils.createRealMatrix(data);
109             RealMatrix h = checkAEqualPHPt(m);
110             checkHessenbergForm(h);
111         }
112     }
113 
114     @Test
115     public void testRandomDataNormalDistribution() {
116         for (int run = 0; run < 100; run++) {
117             Random r = new Random(System.currentTimeMillis());
118             RandomDataGenerator gen = new RandomDataGenerator(100);
119 
120             // matrix size
121             int size = r.nextInt(20) + 4;
122 
123             double[][] data = new double[size][size];
124             for (int i = 0; i < size; i++) {
125                 for (int j = 0; j < size; j++) {
126                     data[i][j] = gen.nextNormal(0.0, r.nextDouble() * 5);
127                 }
128             }
129 
130             RealMatrix m = MatrixUtils.createRealMatrix(data);
131             RealMatrix h = checkAEqualPHPt(m);
132             checkHessenbergForm(h);
133         }
134     }
135 
136     @Test
137     public void testMatricesValues5() {
138         checkMatricesValues(testSquare5,
139                             new double[][] {
140                                 { 1.0,  0.0,                0.0,                0.0,                0.0               },
141                                 { 0.0, -0.182574185835055,  0.784218758628863,  0.395029040913988, -0.442289115981669 },
142                                 { 0.0, -0.365148371670111, -0.337950625265477, -0.374110794088820, -0.782621974707823 },
143                                 { 0.0, -0.547722557505166,  0.402941130124223, -0.626468266309003,  0.381019628053472 },
144                                 { 0.0, -0.730296743340221, -0.329285224617644,  0.558149336547665,  0.216118545309225 }
145                             },
146                             new double[][] {
147                                 {  5.0,              -3.65148371670111,  2.59962019434982, -0.237003414680848, -3.13886458663398  },
148                                 { -5.47722557505166,  6.9,              -2.29164066120599,  0.207283564429169,  0.703858369151728 },
149                                 {  0.0,              -4.21386600008432,  2.30555659846067,  2.74935928725112,   0.857569835914113 },
150                                 {  0.0,               0.0,               2.86406180891882, -1.11582249161595,   0.817995267184158 },
151                                 {  0.0,               0.0,               0.0,               0.683518597386085,  1.91026589315528  }
152                             });
153     }
154 
155     @Test
156     public void testMatricesValues3() {
157         checkMatricesValues(testSquare3,
158                             new double[][] {
159                                 {  1.0,  0.0,               0.0               },
160                                 {  0.0, -0.707106781186547, 0.707106781186547 },
161                                 {  0.0,  0.707106781186547, 0.707106781186548 },
162                             },
163                             new double[][] {
164                                 {  2.0,              1.41421356237309,  0.0 },
165                                 {  1.41421356237310, 2.0,              -1.0 },
166                                 {  0.0,              1.0,               2.0 },
167                             });
168     }
169 
170     ///////////////////////////////////////////////////////////////////////////
171     // Test helpers
172     ///////////////////////////////////////////////////////////////////////////
173 
174     private RealMatrix checkAEqualPHPt(RealMatrix matrix) {
175         HessenbergTransformer transformer = new HessenbergTransformer(matrix);
176         RealMatrix p  = transformer.getP();
177         RealMatrix pT = transformer.getPT();
178         RealMatrix h  = transformer.getH();
179 
180         RealMatrix result = p.multiply(h).multiply(pT);
181         double norm = result.subtract(matrix).getNorm1();
182         Assert.assertEquals(0, norm, 1.0e-10);
183 
184         for (int i = 0; i < matrix.getRowDimension(); ++i) {
185             for (int j = 0; j < matrix.getColumnDimension(); ++j) {
186                 if (i > j + 1) {
187                     Assert.assertEquals(matrix.getEntry(i, j), result.getEntry(i, j), 1.0e-12);
188                 }
189             }
190         }
191 
192         return transformer.getH();
193     }
194 
195     private void checkOrthogonal(RealMatrix m) {
196         RealMatrix mTm = m.transposeMultiply(m);
197         RealMatrix id  = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
198         Assert.assertEquals(0, mTm.subtract(id).getNorm1(), 1.0e-14);
199     }
200 
201     private void checkHessenbergForm(RealMatrix m) {
202         final int rows = m.getRowDimension();
203         final int cols = m.getColumnDimension();
204         for (int i = 0; i < rows; ++i) {
205             for (int j = 0; j < cols; ++j) {
206                 if (i > j + 1) {
207                     Assert.assertEquals(0, m.getEntry(i, j), 1.0e-16);
208                 }
209             }
210         }
211     }
212 
213     private void checkMatricesValues(double[][] matrix, double[][] pRef, double[][] hRef) {
214         HessenbergTransformer transformer =
215             new HessenbergTransformer(MatrixUtils.createRealMatrix(matrix));
216 
217         // check values against known references
218         RealMatrix p = transformer.getP();
219         Assert.assertEquals(0, p.subtract(MatrixUtils.createRealMatrix(pRef)).getNorm1(), 1.0e-14);
220 
221         RealMatrix h = transformer.getH();
222         Assert.assertEquals(0, h.subtract(MatrixUtils.createRealMatrix(hRef)).getNorm1(), 1.3e-14);
223 
224         // check the same cached instance is returned the second time
225         Assert.assertTrue(p == transformer.getP());
226         Assert.assertTrue(h == transformer.getH());
227     }
228 }