source: src/main/java/agents/anac/y2019/harddealer/math3/linear/BiDiagonalTransformer.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: 13.1 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;
21
22
23/**
24 * Class transforming any matrix to bi-diagonal shape.
25 * <p>Any m &times; n matrix A can be written as the product of three matrices:
26 * A = U &times; B &times; V<sup>T</sup> with U an m &times; m orthogonal matrix,
27 * B an m &times; n bi-diagonal matrix (lower diagonal if m &lt; n, upper diagonal
28 * otherwise), and V an n &times; n orthogonal matrix.</p>
29 * <p>Transformation to bi-diagonal shape is often not a goal by itself, but it is
30 * an intermediate step in more general decomposition algorithms like {@link
31 * SingularValueDecomposition Singular Value Decomposition}. This class is therefore
32 * intended for internal use by the library and is not public. As a consequence of
33 * this explicitly limited scope, many methods directly returns references to
34 * internal arrays, not copies.</p>
35 * @since 2.0
36 */
37class BiDiagonalTransformer {
38
39 /** Householder vectors. */
40 private final double householderVectors[][];
41
42 /** Main diagonal. */
43 private final double[] main;
44
45 /** Secondary diagonal. */
46 private final double[] secondary;
47
48 /** Cached value of U. */
49 private RealMatrix cachedU;
50
51 /** Cached value of B. */
52 private RealMatrix cachedB;
53
54 /** Cached value of V. */
55 private RealMatrix cachedV;
56
57 /**
58 * Build the transformation to bi-diagonal shape of a matrix.
59 * @param matrix the matrix to transform.
60 */
61 BiDiagonalTransformer(RealMatrix matrix) {
62
63 final int m = matrix.getRowDimension();
64 final int n = matrix.getColumnDimension();
65 final int p = FastMath.min(m, n);
66 householderVectors = matrix.getData();
67 main = new double[p];
68 secondary = new double[p - 1];
69 cachedU = null;
70 cachedB = null;
71 cachedV = null;
72
73 // transform matrix
74 if (m >= n) {
75 transformToUpperBiDiagonal();
76 } else {
77 transformToLowerBiDiagonal();
78 }
79
80 }
81
82 /**
83 * Returns the matrix U of the transform.
84 * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
85 * @return the U matrix
86 */
87 public RealMatrix getU() {
88
89 if (cachedU == null) {
90
91 final int m = householderVectors.length;
92 final int n = householderVectors[0].length;
93 final int p = main.length;
94 final int diagOffset = (m >= n) ? 0 : 1;
95 final double[] diagonal = (m >= n) ? main : secondary;
96 double[][] ua = new double[m][m];
97
98 // fill up the part of the matrix not affected by Householder transforms
99 for (int k = m - 1; k >= p; --k) {
100 ua[k][k] = 1;
101 }
102
103 // build up first part of the matrix by applying Householder transforms
104 for (int k = p - 1; k >= diagOffset; --k) {
105 final double[] hK = householderVectors[k];
106 ua[k][k] = 1;
107 if (hK[k - diagOffset] != 0.0) {
108 for (int j = k; j < m; ++j) {
109 double alpha = 0;
110 for (int i = k; i < m; ++i) {
111 alpha -= ua[i][j] * householderVectors[i][k - diagOffset];
112 }
113 alpha /= diagonal[k - diagOffset] * hK[k - diagOffset];
114
115 for (int i = k; i < m; ++i) {
116 ua[i][j] += -alpha * householderVectors[i][k - diagOffset];
117 }
118 }
119 }
120 }
121 if (diagOffset > 0) {
122 ua[0][0] = 1;
123 }
124 cachedU = MatrixUtils.createRealMatrix(ua);
125 }
126
127 // return the cached matrix
128 return cachedU;
129
130 }
131
132 /**
133 * Returns the bi-diagonal matrix B of the transform.
134 * @return the B matrix
135 */
136 public RealMatrix getB() {
137
138 if (cachedB == null) {
139
140 final int m = householderVectors.length;
141 final int n = householderVectors[0].length;
142 double[][] ba = new double[m][n];
143 for (int i = 0; i < main.length; ++i) {
144 ba[i][i] = main[i];
145 if (m < n) {
146 if (i > 0) {
147 ba[i][i-1] = secondary[i - 1];
148 }
149 } else {
150 if (i < main.length - 1) {
151 ba[i][i+1] = secondary[i];
152 }
153 }
154 }
155 cachedB = MatrixUtils.createRealMatrix(ba);
156 }
157
158 // return the cached matrix
159 return cachedB;
160
161 }
162
163 /**
164 * Returns the matrix V of the transform.
165 * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
166 * @return the V matrix
167 */
168 public RealMatrix getV() {
169
170 if (cachedV == null) {
171
172 final int m = householderVectors.length;
173 final int n = householderVectors[0].length;
174 final int p = main.length;
175 final int diagOffset = (m >= n) ? 1 : 0;
176 final double[] diagonal = (m >= n) ? secondary : main;
177 double[][] va = new double[n][n];
178
179 // fill up the part of the matrix not affected by Householder transforms
180 for (int k = n - 1; k >= p; --k) {
181 va[k][k] = 1;
182 }
183
184 // build up first part of the matrix by applying Householder transforms
185 for (int k = p - 1; k >= diagOffset; --k) {
186 final double[] hK = householderVectors[k - diagOffset];
187 va[k][k] = 1;
188 if (hK[k] != 0.0) {
189 for (int j = k; j < n; ++j) {
190 double beta = 0;
191 for (int i = k; i < n; ++i) {
192 beta -= va[i][j] * hK[i];
193 }
194 beta /= diagonal[k - diagOffset] * hK[k];
195
196 for (int i = k; i < n; ++i) {
197 va[i][j] += -beta * hK[i];
198 }
199 }
200 }
201 }
202 if (diagOffset > 0) {
203 va[0][0] = 1;
204 }
205 cachedV = MatrixUtils.createRealMatrix(va);
206 }
207
208 // return the cached matrix
209 return cachedV;
210
211 }
212
213 /**
214 * Get the Householder vectors of the transform.
215 * <p>Note that since this class is only intended for internal use,
216 * it returns directly a reference to its internal arrays, not a copy.</p>
217 * @return the main diagonal elements of the B matrix
218 */
219 double[][] getHouseholderVectorsRef() {
220 return householderVectors;
221 }
222
223 /**
224 * Get the main diagonal elements of the matrix B of the transform.
225 * <p>Note that since this class is only intended for internal use,
226 * it returns directly a reference to its internal arrays, not a copy.</p>
227 * @return the main diagonal elements of the B matrix
228 */
229 double[] getMainDiagonalRef() {
230 return main;
231 }
232
233 /**
234 * Get the secondary diagonal elements of the matrix B of the transform.
235 * <p>Note that since this class is only intended for internal use,
236 * it returns directly a reference to its internal arrays, not a copy.</p>
237 * @return the secondary diagonal elements of the B matrix
238 */
239 double[] getSecondaryDiagonalRef() {
240 return secondary;
241 }
242
243 /**
244 * Check if the matrix is transformed to upper bi-diagonal.
245 * @return true if the matrix is transformed to upper bi-diagonal
246 */
247 boolean isUpperBiDiagonal() {
248 return householderVectors.length >= householderVectors[0].length;
249 }
250
251 /**
252 * Transform original matrix to upper bi-diagonal form.
253 * <p>Transformation is done using alternate Householder transforms
254 * on columns and rows.</p>
255 */
256 private void transformToUpperBiDiagonal() {
257
258 final int m = householderVectors.length;
259 final int n = householderVectors[0].length;
260 for (int k = 0; k < n; k++) {
261
262 //zero-out a column
263 double xNormSqr = 0;
264 for (int i = k; i < m; ++i) {
265 final double c = householderVectors[i][k];
266 xNormSqr += c * c;
267 }
268 final double[] hK = householderVectors[k];
269 final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
270 main[k] = a;
271 if (a != 0.0) {
272 hK[k] -= a;
273 for (int j = k + 1; j < n; ++j) {
274 double alpha = 0;
275 for (int i = k; i < m; ++i) {
276 final double[] hI = householderVectors[i];
277 alpha -= hI[j] * hI[k];
278 }
279 alpha /= a * householderVectors[k][k];
280 for (int i = k; i < m; ++i) {
281 final double[] hI = householderVectors[i];
282 hI[j] -= alpha * hI[k];
283 }
284 }
285 }
286
287 if (k < n - 1) {
288 //zero-out a row
289 xNormSqr = 0;
290 for (int j = k + 1; j < n; ++j) {
291 final double c = hK[j];
292 xNormSqr += c * c;
293 }
294 final double b = (hK[k + 1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
295 secondary[k] = b;
296 if (b != 0.0) {
297 hK[k + 1] -= b;
298 for (int i = k + 1; i < m; ++i) {
299 final double[] hI = householderVectors[i];
300 double beta = 0;
301 for (int j = k + 1; j < n; ++j) {
302 beta -= hI[j] * hK[j];
303 }
304 beta /= b * hK[k + 1];
305 for (int j = k + 1; j < n; ++j) {
306 hI[j] -= beta * hK[j];
307 }
308 }
309 }
310 }
311
312 }
313 }
314
315 /**
316 * Transform original matrix to lower bi-diagonal form.
317 * <p>Transformation is done using alternate Householder transforms
318 * on rows and columns.</p>
319 */
320 private void transformToLowerBiDiagonal() {
321
322 final int m = householderVectors.length;
323 final int n = householderVectors[0].length;
324 for (int k = 0; k < m; k++) {
325
326 //zero-out a row
327 final double[] hK = householderVectors[k];
328 double xNormSqr = 0;
329 for (int j = k; j < n; ++j) {
330 final double c = hK[j];
331 xNormSqr += c * c;
332 }
333 final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
334 main[k] = a;
335 if (a != 0.0) {
336 hK[k] -= a;
337 for (int i = k + 1; i < m; ++i) {
338 final double[] hI = householderVectors[i];
339 double alpha = 0;
340 for (int j = k; j < n; ++j) {
341 alpha -= hI[j] * hK[j];
342 }
343 alpha /= a * householderVectors[k][k];
344 for (int j = k; j < n; ++j) {
345 hI[j] -= alpha * hK[j];
346 }
347 }
348 }
349
350 if (k < m - 1) {
351 //zero-out a column
352 final double[] hKp1 = householderVectors[k + 1];
353 xNormSqr = 0;
354 for (int i = k + 1; i < m; ++i) {
355 final double c = householderVectors[i][k];
356 xNormSqr += c * c;
357 }
358 final double b = (hKp1[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
359 secondary[k] = b;
360 if (b != 0.0) {
361 hKp1[k] -= b;
362 for (int j = k + 1; j < n; ++j) {
363 double beta = 0;
364 for (int i = k + 1; i < m; ++i) {
365 final double[] hI = householderVectors[i];
366 beta -= hI[j] * hI[k];
367 }
368 beta /= b * hKp1[k];
369 for (int i = k + 1; i < m; ++i) {
370 final double[] hI = householderVectors[i];
371 hI[j] -= beta * hI[k];
372 }
373 }
374 }
375 }
376
377 }
378 }
379
380}
Note: See TracBrowser for help on using the repository browser.