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