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    *      http://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  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 SchurTransformerTest {
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      // from http://eigen.tuxfamily.org/dox/classEigen_1_1RealSchur.html
45      private double[][] testRandom = {
46              {  0.680, -0.3300, -0.2700, -0.717, -0.687,  0.0259 },
47              { -0.211,  0.5360,  0.0268,  0.214, -0.198,  0.6780 },
48              {  0.566, -0.4440,  0.9040, -0.967, -0.740,  0.2250 },
49              {  0.597,  0.1080,  0.8320, -0.514, -0.782, -0.4080 },
50              {  0.823, -0.0452,  0.2710, -0.726,  0.998,  0.2750 },
51              { -0.605,  0.2580,  0.4350,  0.608, -0.563,  0.0486 }
52      };
53  
54      @Test
55      public void testNonSquare() {
56          try {
57              new SchurTransformer(MatrixUtils.createRealMatrix(new double[3][2]));
58              Assert.fail("an exception should have been thrown");
59          } catch (NonSquareMatrixException ime) {
60              // expected behavior
61          }
62      }
63  
64      @Test
65      public void testAEqualPTPt() {
66          checkAEqualPTPt(MatrixUtils.createRealMatrix(testSquare5));
67          checkAEqualPTPt(MatrixUtils.createRealMatrix(testSquare3));
68          checkAEqualPTPt(MatrixUtils.createRealMatrix(testRandom));
69     }
70  
71      @Test
72      public void testPOrthogonal() {
73          checkOrthogonal(new SchurTransformer(MatrixUtils.createRealMatrix(testSquare5)).getP());
74          checkOrthogonal(new SchurTransformer(MatrixUtils.createRealMatrix(testSquare3)).getP());
75          checkOrthogonal(new SchurTransformer(MatrixUtils.createRealMatrix(testRandom)).getP());
76      }
77  
78      @Test
79      public void testPTOrthogonal() {
80          checkOrthogonal(new SchurTransformer(MatrixUtils.createRealMatrix(testSquare5)).getPT());
81          checkOrthogonal(new SchurTransformer(MatrixUtils.createRealMatrix(testSquare3)).getPT());
82          checkOrthogonal(new SchurTransformer(MatrixUtils.createRealMatrix(testRandom)).getPT());
83      }
84  
85      @Test
86      public void testSchurForm() {
87          checkSchurForm(new SchurTransformer(MatrixUtils.createRealMatrix(testSquare5)).getT());
88          checkSchurForm(new SchurTransformer(MatrixUtils.createRealMatrix(testSquare3)).getT());
89          checkSchurForm(new SchurTransformer(MatrixUtils.createRealMatrix(testRandom)).getT());
90      }
91  
92      @Test
93      public void testRandomData() {
94          for (int run = 0; run < 100; run++) {
95              Random r = new Random(System.currentTimeMillis());
96  
97              // matrix size
98              int size = r.nextInt(20) + 4;
99  
100             double[][] data = new double[size][size];
101             for (int i = 0; i < size; i++) {
102                 for (int j = 0; j < size; j++) {
103                     data[i][j] = r.nextInt(100);
104                 }
105             }
106 
107             RealMatrix m = MatrixUtils.createRealMatrix(data);
108             RealMatrix s = checkAEqualPTPt(m);
109             checkSchurForm(s);
110         }
111     }
112 
113     @Test
114     public void testRandomDataNormalDistribution() {
115         for (int run = 0; run < 100; run++) {
116             Random r = new Random(System.currentTimeMillis());
117             ContinuousDistribution.Sampler dist
118                 = NormalDistribution.of(0.0, r.nextDouble() * 5).createSampler(RandomSource.WELL_512_A.create(64925784252L));
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] = dist.sample();
127                 }
128             }
129 
130             RealMatrix m = MatrixUtils.createRealMatrix(data);
131             RealMatrix s = checkAEqualPTPt(m);
132             checkSchurForm(s);
133         }
134     }
135 
136     @Test
137     public void testMath848() {
138         double[][] data = {
139                 { 0.1849449280, -0.0646971046,  0.0774755812, -0.0969651755, -0.0692648806,  0.3282344352, -0.0177423074,  0.2063136340},
140                 {-0.0742700134, -0.0289063030, -0.0017269460, -0.0375550146, -0.0487737922, -0.2616837868, -0.0821201295, -0.2530000167},
141                 { 0.2549910127,  0.0995733692, -0.0009718388,  0.0149282808,  0.1791878897, -0.0823182816,  0.0582629256,  0.3219545182},
142                 {-0.0694747557, -0.1880649148, -0.2740630911,  0.0720096468, -0.1800836914, -0.3518996425,  0.2486747833,  0.6257938167},
143                 { 0.0536360918, -0.1339297778,  0.2241579764, -0.0195327484, -0.0054103808,  0.0347564518,  0.5120802482, -0.0329902864},
144                 {-0.5933332356, -0.2488721082,  0.2357173629,  0.0177285473,  0.0856630593, -0.3567126300, -0.1600668126, -0.1010899621},
145                 {-0.0514349819, -0.0854319435,  0.1125050061,  0.0063453560, -0.2250000688, -0.2209343090,  0.1964623477, -0.1512329924},
146                 { 0.0197395947, -0.1997170581, -0.1425959019, -0.2749477910, -0.0969467073,  0.0603688520, -0.2826905192,  0.1794315473}};
147         RealMatrix m = MatrixUtils.createRealMatrix(data);
148         RealMatrix s = checkAEqualPTPt(m);
149         checkSchurForm(s);
150     }
151 
152     ///////////////////////////////////////////////////////////////////////////
153     // Test helpers
154     ///////////////////////////////////////////////////////////////////////////
155 
156     private RealMatrix checkAEqualPTPt(RealMatrix matrix) {
157         SchurTransformer transformer = new SchurTransformer(matrix);
158         RealMatrix p  = transformer.getP();
159         RealMatrix t  = transformer.getT();
160         RealMatrix pT = transformer.getPT();
161 
162         RealMatrix result = p.multiply(t).multiply(pT);
163 
164         double norm = result.subtract(matrix).getNorm();
165         Assert.assertEquals(0, norm, 1.0e-9);
166 
167         return t;
168     }
169 
170     private void checkOrthogonal(RealMatrix m) {
171         RealMatrix mTm = m.transpose().multiply(m);
172         RealMatrix id  = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
173         Assert.assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);
174     }
175 
176     private void checkSchurForm(final RealMatrix m) {
177         final int rows = m.getRowDimension();
178         final int cols = m.getColumnDimension();
179         for (int i = 0; i < rows; ++i) {
180             for (int j = 0; j < cols; ++j) {
181                 if (i > j + 1) {
182                     Assert.assertEquals(0, m.getEntry(i, j), 1.0e-16);
183                 }
184             }
185         }
186     }
187 
188     @SuppressWarnings("unused")
189     private void checkMatricesValues(double[][] matrix, double[][] pRef, double[][] hRef) {
190 
191         SchurTransformer transformer =
192             new SchurTransformer(MatrixUtils.createRealMatrix(matrix));
193 
194         // check values against known references
195         RealMatrix p = transformer.getP();
196         Assert.assertEquals(0, p.subtract(MatrixUtils.createRealMatrix(pRef)).getNorm(), 1.0e-14);
197 
198         RealMatrix t = transformer.getT();
199         Assert.assertEquals(0, t.subtract(MatrixUtils.createRealMatrix(hRef)).getNorm(), 1.0e-14);
200 
201         // check the same cached instance is returned the second time
202         Assert.assertSame(p, transformer.getP());
203         Assert.assertSame(t, transformer.getT());
204     }
205 }