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 org.apache.commons.math4.core.jdkmath.JdkMath;
21
22
23 /**
24 * Class transforming any matrix to bi-diagonal shape.
25 * <p>Any m × n matrix A can be written as the product of three matrices:
26 * A = U × B × V<sup>T</sup> with U an m × m orthogonal matrix,
27 * B an m × n bi-diagonal matrix (lower diagonal if m < n, upper diagonal
28 * otherwise), and V an n × 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 */
37 class 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 = JdkMath.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 * Returns the matrix U of the transform.
83 * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
84 * @return the U matrix
85 */
86 public RealMatrix getU() {
87
88 if (cachedU == null) {
89
90 final int m = householderVectors.length;
91 final int n = householderVectors[0].length;
92 final int p = main.length;
93 final int diagOffset = (m >= n) ? 0 : 1;
94 final double[] diagonal = (m >= n) ? main : secondary;
95 double[][] ua = new double[m][m];
96
97 // fill up the part of the matrix not affected by Householder transforms
98 for (int k = m - 1; k >= p; --k) {
99 ua[k][k] = 1;
100 }
101
102 // build up first part of the matrix by applying Householder transforms
103 for (int k = p - 1; k >= diagOffset; --k) {
104 final double[] hK = householderVectors[k];
105 ua[k][k] = 1;
106 if (hK[k - diagOffset] != 0.0) {
107 for (int j = k; j < m; ++j) {
108 double alpha = 0;
109 for (int i = k; i < m; ++i) {
110 alpha -= ua[i][j] * householderVectors[i][k - diagOffset];
111 }
112 alpha /= diagonal[k - diagOffset] * hK[k - diagOffset];
113
114 for (int i = k; i < m; ++i) {
115 ua[i][j] += -alpha * householderVectors[i][k - diagOffset];
116 }
117 }
118 }
119 }
120 if (diagOffset > 0) {
121 ua[0][0] = 1;
122 }
123 cachedU = MatrixUtils.createRealMatrix(ua);
124 }
125
126 // return the cached matrix
127 return cachedU;
128 }
129
130 /**
131 * Returns the bi-diagonal matrix B of the transform.
132 * @return the B matrix
133 */
134 public RealMatrix getB() {
135
136 if (cachedB == null) {
137
138 final int m = householderVectors.length;
139 final int n = householderVectors[0].length;
140 double[][] ba = new double[m][n];
141 for (int i = 0; i < main.length; ++i) {
142 ba[i][i] = main[i];
143 if (m < n) {
144 if (i > 0) {
145 ba[i][i-1] = secondary[i - 1];
146 }
147 } else {
148 if (i < main.length - 1) {
149 ba[i][i+1] = secondary[i];
150 }
151 }
152 }
153 cachedB = MatrixUtils.createRealMatrix(ba);
154 }
155
156 // return the cached matrix
157 return cachedB;
158 }
159
160 /**
161 * Returns the matrix V of the transform.
162 * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
163 * @return the V matrix
164 */
165 public RealMatrix getV() {
166
167 if (cachedV == null) {
168
169 final int m = householderVectors.length;
170 final int n = householderVectors[0].length;
171 final int p = main.length;
172 final int diagOffset = (m >= n) ? 1 : 0;
173 final double[] diagonal = (m >= n) ? secondary : main;
174 double[][] va = new double[n][n];
175
176 // fill up the part of the matrix not affected by Householder transforms
177 for (int k = n - 1; k >= p; --k) {
178 va[k][k] = 1;
179 }
180
181 // build up first part of the matrix by applying Householder transforms
182 for (int k = p - 1; k >= diagOffset; --k) {
183 final double[] hK = householderVectors[k - diagOffset];
184 va[k][k] = 1;
185 if (hK[k] != 0.0) {
186 for (int j = k; j < n; ++j) {
187 double beta = 0;
188 for (int i = k; i < n; ++i) {
189 beta -= va[i][j] * hK[i];
190 }
191 beta /= diagonal[k - diagOffset] * hK[k];
192
193 for (int i = k; i < n; ++i) {
194 va[i][j] += -beta * hK[i];
195 }
196 }
197 }
198 }
199 if (diagOffset > 0) {
200 va[0][0] = 1;
201 }
202 cachedV = MatrixUtils.createRealMatrix(va);
203 }
204
205 // return the cached matrix
206 return cachedV;
207 }
208
209 /**
210 * Get the Householder vectors of the transform.
211 * <p>Note that since this class is only intended for internal use,
212 * it returns directly a reference to its internal arrays, not a copy.</p>
213 * @return the main diagonal elements of the B matrix
214 */
215 double[][] getHouseholderVectorsRef() {
216 return householderVectors;
217 }
218
219 /**
220 * Get the main diagonal elements of the matrix B of the transform.
221 * <p>Note that since this class is only intended for internal use,
222 * it returns directly a reference to its internal arrays, not a copy.</p>
223 * @return the main diagonal elements of the B matrix
224 */
225 double[] getMainDiagonalRef() {
226 return main;
227 }
228
229 /**
230 * Get the secondary diagonal elements of the matrix B of the transform.
231 * <p>Note that since this class is only intended for internal use,
232 * it returns directly a reference to its internal arrays, not a copy.</p>
233 * @return the secondary diagonal elements of the B matrix
234 */
235 double[] getSecondaryDiagonalRef() {
236 return secondary;
237 }
238
239 /**
240 * Check if the matrix is transformed to upper bi-diagonal.
241 * @return true if the matrix is transformed to upper bi-diagonal
242 */
243 boolean isUpperBiDiagonal() {
244 return householderVectors.length >= householderVectors[0].length;
245 }
246
247 /**
248 * Transform original matrix to upper bi-diagonal form.
249 * <p>Transformation is done using alternate Householder transforms
250 * on columns and rows.</p>
251 */
252 private void transformToUpperBiDiagonal() {
253
254 final int m = householderVectors.length;
255 final int n = householderVectors[0].length;
256 for (int k = 0; k < n; k++) {
257
258 //zero-out a column
259 double xNormSqr = 0;
260 for (int i = k; i < m; ++i) {
261 final double c = householderVectors[i][k];
262 xNormSqr += c * c;
263 }
264 final double[] hK = householderVectors[k];
265 final double a = (hK[k] > 0) ? -JdkMath.sqrt(xNormSqr) : JdkMath.sqrt(xNormSqr);
266 main[k] = a;
267 if (a != 0.0) {
268 hK[k] -= a;
269 for (int j = k + 1; j < n; ++j) {
270 double alpha = 0;
271 for (int i = k; i < m; ++i) {
272 final double[] hI = householderVectors[i];
273 alpha -= hI[j] * hI[k];
274 }
275 alpha /= a * householderVectors[k][k];
276 for (int i = k; i < m; ++i) {
277 final double[] hI = householderVectors[i];
278 hI[j] -= alpha * hI[k];
279 }
280 }
281 }
282
283 if (k < n - 1) {
284 //zero-out a row
285 xNormSqr = 0;
286 for (int j = k + 1; j < n; ++j) {
287 final double c = hK[j];
288 xNormSqr += c * c;
289 }
290 final double b = (hK[k + 1] > 0) ? -JdkMath.sqrt(xNormSqr) : JdkMath.sqrt(xNormSqr);
291 secondary[k] = b;
292 if (b != 0.0) {
293 hK[k + 1] -= b;
294 for (int i = k + 1; i < m; ++i) {
295 final double[] hI = householderVectors[i];
296 double beta = 0;
297 for (int j = k + 1; j < n; ++j) {
298 beta -= hI[j] * hK[j];
299 }
300 beta /= b * hK[k + 1];
301 for (int j = k + 1; j < n; ++j) {
302 hI[j] -= beta * hK[j];
303 }
304 }
305 }
306 }
307 }
308 }
309
310 /**
311 * Transform original matrix to lower bi-diagonal form.
312 * <p>Transformation is done using alternate Householder transforms
313 * on rows and columns.</p>
314 */
315 private void transformToLowerBiDiagonal() {
316
317 final int m = householderVectors.length;
318 final int n = householderVectors[0].length;
319 for (int k = 0; k < m; k++) {
320
321 //zero-out a row
322 final double[] hK = householderVectors[k];
323 double xNormSqr = 0;
324 for (int j = k; j < n; ++j) {
325 final double c = hK[j];
326 xNormSqr += c * c;
327 }
328 final double a = (hK[k] > 0) ? -JdkMath.sqrt(xNormSqr) : JdkMath.sqrt(xNormSqr);
329 main[k] = a;
330 if (a != 0.0) {
331 hK[k] -= a;
332 for (int i = k + 1; i < m; ++i) {
333 final double[] hI = householderVectors[i];
334 double alpha = 0;
335 for (int j = k; j < n; ++j) {
336 alpha -= hI[j] * hK[j];
337 }
338 alpha /= a * householderVectors[k][k];
339 for (int j = k; j < n; ++j) {
340 hI[j] -= alpha * hK[j];
341 }
342 }
343 }
344
345 if (k < m - 1) {
346 //zero-out a column
347 final double[] hKp1 = householderVectors[k + 1];
348 xNormSqr = 0;
349 for (int i = k + 1; i < m; ++i) {
350 final double c = householderVectors[i][k];
351 xNormSqr += c * c;
352 }
353 final double b = (hKp1[k] > 0) ? -JdkMath.sqrt(xNormSqr) : JdkMath.sqrt(xNormSqr);
354 secondary[k] = b;
355 if (b != 0.0) {
356 hKp1[k] -= b;
357 for (int j = k + 1; j < n; ++j) {
358 double beta = 0;
359 for (int i = k + 1; i < m; ++i) {
360 final double[] hI = householderVectors[i];
361 beta -= hI[j] * hI[k];
362 }
363 beta /= b * hKp1[k];
364 for (int i = k + 1; i < m; ++i) {
365 final double[] hI = householderVectors[i];
366 hI[j] -= beta * hI[k];
367 }
368 }
369 }
370 }
371 }
372 }
373 }