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