source: src/main/java/agents/anac/y2019/harddealer/math3/linear/HessenbergTransformer.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 8.0 KB
Line 
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
18package agents.anac.y2019.harddealer.math3.linear;
19
20import agents.anac.y2019.harddealer.math3.util.FastMath;
21import agents.anac.y2019.harddealer.math3.util.Precision;
22
23/**
24 * Class transforming a general real matrix to Hessenberg form.
25 * <p>A m &times; m matrix A can be written as the product of three matrices: A = P
26 * &times; H &times; P<sup>T</sup> with P an orthogonal matrix and H a Hessenberg
27 * matrix. Both P and H are m &times; m matrices.</p>
28 * <p>Transformation to Hessenberg form is often not a goal by itself, but it is an
29 * intermediate step in more general decomposition algorithms like
30 * {@link EigenDecomposition eigen decomposition}. This class is therefore
31 * intended for internal use by the library and is not public. As a consequence
32 * of this explicitly limited scope, many methods directly returns references to
33 * internal arrays, not copies.</p>
34 * <p>This class is based on the method orthes in class EigenvalueDecomposition
35 * from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
36 *
37 * @see <a href="http://mathworld.wolfram.com/HessenbergDecomposition.html">MathWorld</a>
38 * @see <a href="http://en.wikipedia.org/wiki/Householder_transformation">Householder Transformations</a>
39 * @since 3.1
40 */
41class HessenbergTransformer {
42 /** Householder vectors. */
43 private final double householderVectors[][];
44 /** Temporary storage vector. */
45 private final double ort[];
46 /** Cached value of P. */
47 private RealMatrix cachedP;
48 /** Cached value of Pt. */
49 private RealMatrix cachedPt;
50 /** Cached value of H. */
51 private RealMatrix cachedH;
52
53 /**
54 * Build the transformation to Hessenberg form of a general matrix.
55 *
56 * @param matrix matrix to transform
57 * @throws NonSquareMatrixException if the matrix is not square
58 */
59 HessenbergTransformer(final RealMatrix matrix) {
60 if (!matrix.isSquare()) {
61 throw new NonSquareMatrixException(matrix.getRowDimension(),
62 matrix.getColumnDimension());
63 }
64
65 final int m = matrix.getRowDimension();
66 householderVectors = matrix.getData();
67 ort = new double[m];
68 cachedP = null;
69 cachedPt = null;
70 cachedH = null;
71
72 // transform matrix
73 transform();
74 }
75
76 /**
77 * Returns the matrix P of the transform.
78 * <p>P is an orthogonal matrix, i.e. its inverse is also its transpose.</p>
79 *
80 * @return the P matrix
81 */
82 public RealMatrix getP() {
83 if (cachedP == null) {
84 final int n = householderVectors.length;
85 final int high = n - 1;
86 final double[][] pa = new double[n][n];
87
88 for (int i = 0; i < n; i++) {
89 for (int j = 0; j < n; j++) {
90 pa[i][j] = (i == j) ? 1 : 0;
91 }
92 }
93
94 for (int m = high - 1; m >= 1; m--) {
95 if (householderVectors[m][m - 1] != 0.0) {
96 for (int i = m + 1; i <= high; i++) {
97 ort[i] = householderVectors[i][m - 1];
98 }
99
100 for (int j = m; j <= high; j++) {
101 double g = 0.0;
102
103 for (int i = m; i <= high; i++) {
104 g += ort[i] * pa[i][j];
105 }
106
107 // Double division avoids possible underflow
108 g = (g / ort[m]) / householderVectors[m][m - 1];
109
110 for (int i = m; i <= high; i++) {
111 pa[i][j] += g * ort[i];
112 }
113 }
114 }
115 }
116
117 cachedP = MatrixUtils.createRealMatrix(pa);
118 }
119 return cachedP;
120 }
121
122 /**
123 * Returns the transpose of the matrix P of the transform.
124 * <p>P is an orthogonal matrix, i.e. its inverse is also its transpose.</p>
125 *
126 * @return the transpose of the P matrix
127 */
128 public RealMatrix getPT() {
129 if (cachedPt == null) {
130 cachedPt = getP().transpose();
131 }
132
133 // return the cached matrix
134 return cachedPt;
135 }
136
137 /**
138 * Returns the Hessenberg matrix H of the transform.
139 *
140 * @return the H matrix
141 */
142 public RealMatrix getH() {
143 if (cachedH == null) {
144 final int m = householderVectors.length;
145 final double[][] h = new double[m][m];
146 for (int i = 0; i < m; ++i) {
147 if (i > 0) {
148 // copy the entry of the lower sub-diagonal
149 h[i][i - 1] = householderVectors[i][i - 1];
150 }
151
152 // copy upper triangular part of the matrix
153 for (int j = i; j < m; ++j) {
154 h[i][j] = householderVectors[i][j];
155 }
156 }
157 cachedH = MatrixUtils.createRealMatrix(h);
158 }
159
160 // return the cached matrix
161 return cachedH;
162 }
163
164 /**
165 * Get the Householder vectors of the transform.
166 * <p>Note that since this class is only intended for internal use, it returns
167 * directly a reference to its internal arrays, not a copy.</p>
168 *
169 * @return the main diagonal elements of the B matrix
170 */
171 double[][] getHouseholderVectorsRef() {
172 return householderVectors;
173 }
174
175 /**
176 * Transform original matrix to Hessenberg form.
177 * <p>Transformation is done using Householder transforms.</p>
178 */
179 private void transform() {
180 final int n = householderVectors.length;
181 final int high = n - 1;
182
183 for (int m = 1; m <= high - 1; m++) {
184 // Scale column.
185 double scale = 0;
186 for (int i = m; i <= high; i++) {
187 scale += FastMath.abs(householderVectors[i][m - 1]);
188 }
189
190 if (!Precision.equals(scale, 0)) {
191 // Compute Householder transformation.
192 double h = 0;
193 for (int i = high; i >= m; i--) {
194 ort[i] = householderVectors[i][m - 1] / scale;
195 h += ort[i] * ort[i];
196 }
197 final double g = (ort[m] > 0) ? -FastMath.sqrt(h) : FastMath.sqrt(h);
198
199 h -= ort[m] * g;
200 ort[m] -= g;
201
202 // Apply Householder similarity transformation
203 // H = (I - u*u' / h) * H * (I - u*u' / h)
204
205 for (int j = m; j < n; j++) {
206 double f = 0;
207 for (int i = high; i >= m; i--) {
208 f += ort[i] * householderVectors[i][j];
209 }
210 f /= h;
211 for (int i = m; i <= high; i++) {
212 householderVectors[i][j] -= f * ort[i];
213 }
214 }
215
216 for (int i = 0; i <= high; i++) {
217 double f = 0;
218 for (int j = high; j >= m; j--) {
219 f += ort[j] * householderVectors[i][j];
220 }
221 f /= h;
222 for (int j = m; j <= high; j++) {
223 householderVectors[i][j] -= f * ort[j];
224 }
225 }
226
227 ort[m] = scale * ort[m];
228 householderVectors[m][m - 1] = scale * g;
229 }
230 }
231 }
232}
Note: See TracBrowser for help on using the repository browser.