1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.linear;
19
20 import java.util.Random;
21
22 import org.apache.commons.statistics.distribution.ContinuousDistribution;
23 import org.apache.commons.statistics.distribution.NormalDistribution;
24 import org.apache.commons.rng.simple.RandomSource;
25 import org.junit.Test;
26 import org.junit.Assert;
27
28 public class HessenbergTransformerTest {
29
30 private double[][] testSquare5 = {
31 { 5, 4, 3, 2, 1 },
32 { 1, 4, 0, 3, 3 },
33 { 2, 0, 3, 0, 0 },
34 { 3, 2, 1, 2, 5 },
35 { 4, 2, 1, 4, 1 }
36 };
37
38 private double[][] testSquare3 = {
39 { 2, -1, 1 },
40 { -1, 2, 1 },
41 { 1, -1, 2 }
42 };
43
44
45
46 private double[][] testRandom = {
47 { 0.680, 0.823, -0.4440, -0.2700 },
48 { -0.211, -0.605, 0.1080, 0.0268 },
49 { 0.566, -0.330, -0.0452, 0.9040 },
50 { 0.597, 0.536, 0.2580, 0.8320 }
51 };
52
53 @Test
54 public void testNonSquare() {
55 try {
56 new HessenbergTransformer(MatrixUtils.createRealMatrix(new double[3][2]));
57 Assert.fail("an exception should have been thrown");
58 } catch (NonSquareMatrixException ime) {
59
60 }
61 }
62
63 @Test
64 public void testAEqualPHPt() {
65 checkAEqualPHPt(MatrixUtils.createRealMatrix(testSquare5));
66 checkAEqualPHPt(MatrixUtils.createRealMatrix(testSquare3));
67 checkAEqualPHPt(MatrixUtils.createRealMatrix(testRandom));
68 }
69
70 @Test
71 public void testPOrthogonal() {
72 checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getP());
73 checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getP());
74 }
75
76 @Test
77 public void testPTOrthogonal() {
78 checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getPT());
79 checkOrthogonal(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getPT());
80 }
81
82 @Test
83 public void testHessenbergForm() {
84 checkHessenbergForm(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare5)).getH());
85 checkHessenbergForm(new HessenbergTransformer(MatrixUtils.createRealMatrix(testSquare3)).getH());
86 }
87
88 @Test
89 public void testRandomData() {
90 for (int run = 0; run < 100; run++) {
91 Random r = new Random(System.currentTimeMillis());
92
93
94 int size = r.nextInt(20) + 4;
95
96 double[][] data = new double[size][size];
97 for (int i = 0; i < size; i++) {
98 for (int j = 0; j < size; j++) {
99 data[i][j] = r.nextInt(100);
100 }
101 }
102
103 RealMatrix m = MatrixUtils.createRealMatrix(data);
104 RealMatrix h = checkAEqualPHPt(m);
105 checkHessenbergForm(h);
106 }
107 }
108
109 @Test
110 public void testRandomDataNormalDistribution() {
111 for (int run = 0; run < 100; run++) {
112 Random r = new Random(System.currentTimeMillis());
113 ContinuousDistribution.Sampler dist
114 = NormalDistribution.of(0.0, r.nextDouble() * 5).createSampler(RandomSource.WELL_512_A.create(64925784252L));
115
116
117 int size = r.nextInt(20) + 4;
118
119 double[][] data = new double[size][size];
120 for (int i = 0; i < size; i++) {
121 for (int j = 0; j < size; j++) {
122 data[i][j] = dist.sample();
123 }
124 }
125
126 RealMatrix m = MatrixUtils.createRealMatrix(data);
127 RealMatrix h = checkAEqualPHPt(m);
128 checkHessenbergForm(h);
129 }
130 }
131
132 @Test
133 public void testMatricesValues5() {
134 checkMatricesValues(testSquare5,
135 new double[][] {
136 { 1.0, 0.0, 0.0, 0.0, 0.0 },
137 { 0.0, -0.182574185835055, 0.784218758628863, 0.395029040913988, -0.442289115981669 },
138 { 0.0, -0.365148371670111, -0.337950625265477, -0.374110794088820, -0.782621974707823 },
139 { 0.0, -0.547722557505166, 0.402941130124223, -0.626468266309003, 0.381019628053472 },
140 { 0.0, -0.730296743340221, -0.329285224617644, 0.558149336547665, 0.216118545309225 }
141 },
142 new double[][] {
143 { 5.0, -3.65148371670111, 2.59962019434982, -0.237003414680848, -3.13886458663398 },
144 { -5.47722557505166, 6.9, -2.29164066120599, 0.207283564429169, 0.703858369151728 },
145 { 0.0, -4.21386600008432, 2.30555659846067, 2.74935928725112, 0.857569835914113 },
146 { 0.0, 0.0, 2.86406180891882, -1.11582249161595, 0.817995267184158 },
147 { 0.0, 0.0, 0.0, 0.683518597386085, 1.91026589315528 }
148 });
149 }
150
151 @Test
152 public void testMatricesValues3() {
153 checkMatricesValues(testSquare3,
154 new double[][] {
155 { 1.0, 0.0, 0.0 },
156 { 0.0, -0.707106781186547, 0.707106781186547 },
157 { 0.0, 0.707106781186547, 0.707106781186548 },
158 },
159 new double[][] {
160 { 2.0, 1.41421356237309, 0.0 },
161 { 1.41421356237310, 2.0, -1.0 },
162 { 0.0, 1.0, 2.0 },
163 });
164 }
165
166
167
168
169
170 private RealMatrix checkAEqualPHPt(RealMatrix matrix) {
171 HessenbergTransformer transformer = new HessenbergTransformer(matrix);
172 RealMatrix p = transformer.getP();
173 RealMatrix pT = transformer.getPT();
174 RealMatrix h = transformer.getH();
175
176 RealMatrix result = p.multiply(h).multiply(pT);
177 double norm = result.subtract(matrix).getNorm();
178 Assert.assertEquals(0, norm, 1.0e-10);
179
180 for (int i = 0; i < matrix.getRowDimension(); ++i) {
181 for (int j = 0; j < matrix.getColumnDimension(); ++j) {
182 if (i > j + 1) {
183 Assert.assertEquals(matrix.getEntry(i, j), result.getEntry(i, j), 1.0e-12);
184 }
185 }
186 }
187
188 return transformer.getH();
189 }
190
191 private void checkOrthogonal(RealMatrix m) {
192 RealMatrix mTm = m.transpose().multiply(m);
193 RealMatrix id = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
194 Assert.assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);
195 }
196
197 private void checkHessenbergForm(RealMatrix m) {
198 final int rows = m.getRowDimension();
199 final int cols = m.getColumnDimension();
200 for (int i = 0; i < rows; ++i) {
201 for (int j = 0; j < cols; ++j) {
202 if (i > j + 1) {
203 Assert.assertEquals(0, m.getEntry(i, j), 1.0e-16);
204 }
205 }
206 }
207 }
208
209 private void checkMatricesValues(double[][] matrix, double[][] pRef, double[][] hRef) {
210 HessenbergTransformer transformer =
211 new HessenbergTransformer(MatrixUtils.createRealMatrix(matrix));
212
213
214 RealMatrix p = transformer.getP();
215 Assert.assertEquals(0, p.subtract(MatrixUtils.createRealMatrix(pRef)).getNorm(), 1.0e-14);
216
217 RealMatrix h = transformer.getH();
218 Assert.assertEquals(0, h.subtract(MatrixUtils.createRealMatrix(hRef)).getNorm(), 1.0e-14);
219
220
221 Assert.assertSame(p, transformer.getP());
222 Assert.assertSame(h, transformer.getH());
223 }
224 }