1 /*
2 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.algebra;
17
18 /**
19 * This decomposer computes QR decomposition, which consists on factoring
20 * provided input matrix into an orthogonal matrix (Q) and an upper triangular
21 * matrix (R). In other words, if input matrix is A, then:
22 * A = Q * R
23 */
24 @SuppressWarnings("DuplicatedCode")
25 public class QRDecomposer extends Decomposer {
26
27 /**
28 * Constant defining default round error when determining full rank of
29 * matrices. This value is zero by default.
30 */
31 public static final double DEFAULT_ROUND_ERROR = 1e-8;
32
33 /**
34 * Constant defining minimum allowed round error value when determining full
35 * rank of matrices.
36 */
37 public static final double MIN_ROUND_ERROR = 0.0;
38
39 /**
40 * Internal matrix containing Q factor.
41 */
42 private Matrix q;
43
44 /**
45 * Internal matrix containing R factor.
46 */
47 private Matrix r;
48
49 /**
50 * Boolean indicating whether decomposed matrix is singular.
51 */
52 boolean sing;
53
54 /**
55 * Constructor of this class.
56 */
57 public QRDecomposer() {
58 super();
59 q = r = null;
60 sing = false;
61 }
62
63 /**
64 * Constructor of this class.
65 *
66 * @param inputMatrix Reference to input matrix to be decomposed.
67 */
68 public QRDecomposer(final Matrix inputMatrix) {
69 super(inputMatrix);
70 q = r = null;
71 sing = false;
72 }
73
74 /**
75 * Returns decomposer type corresponding to QR decomposition.
76 *
77 * @return Decomposer type.
78 */
79 @Override
80 public DecomposerType getDecomposerType() {
81 return DecomposerType.QR_DECOMPOSITION;
82 }
83
84 /**
85 * Sets reference to input matrix to be decomposed.
86 *
87 * @param inputMatrix Reference to input matrix to be decomposed.
88 * @throws LockedException Exception thrown if attempting to call this
89 * method while this instance remains locked.
90 */
91 @Override
92 public void setInputMatrix(final Matrix inputMatrix) throws LockedException {
93 super.setInputMatrix(inputMatrix);
94 q = r = null;
95 sing = false;
96 }
97
98 /**
99 * Returns boolean indicating whether decomposition has been computed and
100 * results can be retrieved.
101 * Attempting to retrieve decomposition results when not available, will
102 * probably raise a NotAvailableException.
103 *
104 * @return Boolean indicating whether decomposition has been computed and
105 * results can be retrieved.
106 */
107 @Override
108 public boolean isDecompositionAvailable() {
109 return q != null || r != null;
110 }
111
112 /**
113 * This method computes LU matrix decomposition, which consists on
114 * retrieving two triangular matrices (Lower triangular and Upper
115 * triangular) as a decomposition of provided input matrix.
116 * In other words, if input matrix is A, then: A = L * U
117 * Note: During execution of this method, this instance will be locked,
118 * and hence attempting to set some parameters might raise a
119 * LockedException.
120 * Note: After execution of this method, LU decomposition will be
121 * available and operations such as retrieving L and U matrices or
122 * computing determinants among others will be able to be done.
123 * Attempting to call any of such operations before calling this method
124 * will raise a NotAvailableException because they require computation of
125 * LU decomposition first.
126 *
127 * @throws NotReadyException Exception thrown if attempting to call this
128 * method when this instance is not ready (i.e. no input matrix has been
129 * provided).
130 * @throws LockedException Exception thrown if attempting to call this
131 * method when this instance is not ready (i.e. no input matrix has been
132 * provided).
133 * @throws DecomposerException Exception thrown if for any reason
134 * decomposition fails while executing, like when convergence of results
135 * can not be obtained, etc.
136 */
137 @Override
138 public void decompose() throws NotReadyException, LockedException, DecomposerException {
139
140 if (!isReady()) {
141 throw new NotReadyException();
142 }
143 if (isLocked()) {
144 throw new LockedException();
145 }
146
147 locked = true;
148
149 final var rows = inputMatrix.getRows();
150 final var columns = inputMatrix.getColumns();
151 if (rows < columns) {
152 throw new DecomposerException();
153 }
154
155 double norm;
156 double prod;
157 try {
158 // initialize Q with some random values (within range of 1 to keep
159 // high accuracy
160 q = Matrix.createWithUniformRandomValues(rows, rows, 0.5, 1.0);
161 r = new Matrix(rows, columns);
162 // initialize factor R to zero (because it will be diagonal)
163 r.initialize(0.0);
164 } catch (final WrongSizeException ignore) {
165 // never happens
166 }
167
168 // Copy contents of input matrix to Q matrix (which is bigger than input
169 // matrix) on top-left area.
170 q.setSubmatrix(0, 0, rows - 1, columns - 1, inputMatrix);
171
172 // Construct QR decomposition by means of Gram-Schmidt
173 for (var j = 0; j < rows; j++) {
174 // Find orthogonal base by means of Gram-Schmidt of all rows of Q
175 // Previous columns of Q will contain already normalized vectors,
176 // while column J must be made orthogonal and then normalized
177 for (var k = 0; k < j; k++) {
178 // compute scalar product of previous k-th orthonormal Q column
179 // with current input matrix j-th column
180 prod = 0.0;
181 for (var i = 0; i < rows; i++) {
182 prod += q.getElementAt(i, k) * q.getElementAt(i, j);
183 }
184
185 // Update R factor with obtained dot product at location (k, j)
186 // (only within r limits when j < columns)
187 if (j < columns) {
188 r.setElementAt(k, j, prod);
189 }
190
191 // Update j-th column of Q by subtracting obtained dot product
192 // respect to previous k-th orthonormal column, on the direction
193 // of this latter column.
194 for (var i = 0; i < rows; i++) {
195 q.setElementAt(i, j, q.getElementAt(i, j) -
196 prod * q.getElementAt(i, k));
197 }
198 }
199 // Normalize column j of Q after computing orthogonal Q column to
200 // make it orthonormal
201 norm = 0.0;
202 for (var i = 0; i < rows; i++) {
203 norm += Math.pow(q.getElementAt(i, j), 2.0); // compute norm
204 }
205 norm = Math.sqrt(norm);
206 for (var i = 0; i < rows; i++) {
207 q.setElementAt(i, j, q.getElementAt(i, j) / norm); // normalize
208 }
209
210 // update R factor diagonal with obtained norm (only within r limits
211 // when j < columns)
212 if (j < columns) {
213 r.setElementAt(j, j, norm);
214 }
215 }
216
217 locked = false;
218 }
219
220 /**
221 * Returns boolean indicating whether provided input matrix has full rank
222 * or not.
223 * Squared matrices having full rank also have determinant different from
224 * zero.
225 * Note: Because of rounding errors testing whether a matrix has full rank
226 * or not might obtain unreliable results especially for non-squared
227 * matrices. In such cases, matrices usually tend to be considered as full
228 * rank even when they are not.
229 * Note: DEFAULT_ROUND_ERROR is used to determine whether matrix is full
230 * rank or not
231 *
232 * @return Boolean indicating whether provided input matrix has full rank or
233 * not.
234 * @throws NotAvailableException Exception thrown if attempting to call this
235 * method before computing QR decomposition. To avoid this exception call
236 * decompose() method first.
237 * @throws IllegalArgumentException Exception thrown if provided rounding
238 * error is lower than minimum allowed value (MIN_ROUND_ERROR)
239 * @see #decompose()
240 */
241 public boolean isFullRank() throws NotAvailableException {
242 return isFullRank(DEFAULT_ROUND_ERROR);
243 }
244
245 /**
246 * Returns boolean indicating whether provided input matrix has full rank
247 * or not.
248 * Squared matrices having full rank also have determinant different from
249 * zero.
250 * Note: Because of rounding errors testing whether a matrix has full rank
251 * or not might obtain unreliable results especially for non-squared
252 * matrices. In such cases, matrices usually tend to be considered as full
253 * rank even when they are not.
254 *
255 * @param roundingError Determines the amount of margin given to determine
256 * whether a matrix has full rank or not due to rounding errors. If not
257 * provided, by default rounding error is set to zero, but this value can
258 * be relaxed if needed.
259 * @return Boolean indicating whether provided input matrix has full rank or
260 * not.
261 * @throws NotAvailableException Exception thrown if attempting to call this
262 * method before computing QR decomposition. To avoid this exception call
263 * decompose() method first.
264 * @throws IllegalArgumentException Exception thrown if provided rounding
265 * error is lower than minimum allowed value (MIN_ROUND_ERROR)
266 * @see #decompose()
267 */
268 public boolean isFullRank(final double roundingError) throws NotAvailableException {
269
270 if (!isDecompositionAvailable()) {
271 throw new NotAvailableException();
272 }
273 if (roundingError < MIN_ROUND_ERROR) {
274 throw new IllegalArgumentException();
275 }
276
277 final var rows = inputMatrix.getRows();
278 final var columns = inputMatrix.getColumns();
279 final var minSize = Math.min(rows, columns);
280
281 for (var j = 0; j < minSize; j++) {
282 if (Math.abs(r.getElementAt(j, j)) <= roundingError) {
283 return false;
284 }
285 }
286 return true;
287 }
288
289
290 /**
291 * Returns upper triangular factor matrix.
292 * QR decomposition decomposes input matrix into Q (orthogonal matrix) and
293 * R, which is an upper triangular matrix.
294 *
295 * @return Upper triangular factor matrix.
296 * @throws NotAvailableException Exception thrown if attempting to call this
297 * method before computing QR decomposition. To avoid this exception call
298 * decompose() method first.
299 * @see #decompose()
300 */
301 public Matrix getR() throws NotAvailableException {
302 if (!isDecompositionAvailable()) {
303 throw new NotAvailableException();
304 }
305
306 return r;
307 }
308
309 /**
310 * Returns the economy-sized orthogonal factor matrix.
311 * QR decomposition decomposes input matrix into Q, which is an orthogonal
312 * matrix and R (upper triangular matrix).
313 *
314 * @return Orthogonal factor matrix.
315 * @throws NotAvailableException Exception thrown if attempting to call this
316 * method before computing QR decomposition. To avoid this exception call
317 * decompose() method first.
318 * @see #decompose()
319 */
320 public Matrix getQ() throws NotAvailableException {
321 if (!isDecompositionAvailable()) {
322 throw new NotAvailableException();
323 }
324
325 return q;
326 }
327
328 /**
329 * Solves a linear system of equations of the following form: A * X = B.
330 * Where A is the input matrix provided for QR decomposition, X is the
331 * solution to the system of equations, and B is the parameters vector/
332 * matrix.
333 * Note: This method can be reused for different b vectors/matrices without
334 * having to recompute QR decomposition on the same input matrix.
335 * Note: Provided b matrix must have the same number of rows as provided
336 * input matrix A, otherwise a WrongSizeException will be raised.
337 * Note: Provided input matrix "A" must have at least as many rows as columns,
338 * otherwise a WrongSizeException will be raised as well. For input matrices
339 * having a higher number of rows than columns, the system of equations will
340 * be overdetermined and the least squares solution will be found.
341 * Note: If provided input matrix A is rank deficient, a
342 * RankDeficientMatrixException will be thrown.
343 * Note: In order to execute this method, a QR decomposition must be
344 * available, otherwise a NotAvailableException will be raised. In order to
345 * avoid this exception call decompose() method first.
346 * Note: To solve the linear system of equations DEFAULT_ROUND_ERROR is used
347 * Note: Solution of linear system of equations is stored in result matrix,
348 * and result matrix is resized if needed.
349 *
350 * @param b Parameters matrix that determines a linear system of equations.
351 * Provided matrix must have the same number of rows as provided input
352 * matrix for QR decomposition. Besides, each column on parameters matrix
353 * will represent a new system of equations, whose solution will be returned
354 * on appropriate column as an output of this method.
355 * @param result Matrix containing solution of linear system of equations on
356 * each column for each column of provided parameters matrix b.
357 * @throws NotAvailableException Exception thrown if attempting to call this
358 * method before computing QR decomposition. To avoid this exception call
359 * decompose() method first.
360 * @throws WrongSizeException Exception thrown if attempting to call this
361 * method using an input matrix with less rows than columns; or if provided
362 * parameters matrix (b) does not have the same number of rows as input
363 * matrix being QR decomposed.
364 * @throws RankDeficientMatrixException Exception thrown if provided input
365 * matrix to be QR decomposed is rank deficient. In this case linear system
366 * of equations cannot be solved.
367 * @see #decompose()
368 */
369 public void solve(final Matrix b, final Matrix result) throws NotAvailableException, WrongSizeException,
370 RankDeficientMatrixException {
371 solve(b, DEFAULT_ROUND_ERROR, result);
372 }
373
374 /**
375 * Solves a linear system of equations of the following form: A * X = B.
376 * Where A is the input matrix provided for QR decomposition, X is the
377 * solution to the system of equations, and B is the parameters vector/
378 * matrix.
379 * Note: This method can be reused for different b vectors/matrices without
380 * having to recompute QR decomposition on the same input matrix.
381 * Note: Provided b matrix must have the same number of rows as provided
382 * input matrix A, otherwise a WrongSizeException will be raised.
383 * Note: Provided input matrix "A" must have at least as many rows as columns,
384 * otherwise a WrongSizeException will be raised as well. For input matrices
385 * having a higher number of rows than columns, the system of equations will
386 * be overdetermined and the least squares solution will be found.
387 * Note: If provided input matrix A is rank deficient, a
388 * RankDeficientMatrixException will be thrown.
389 * Note: In order to execute this method, a QR decomposition must be
390 * available, otherwise a NotAvailableException will be raised. In order to
391 * avoid this exception call decompose() method first.
392 * Note: Solution of linear system of equations is stored in result matrix,
393 * and result matrix is resized if needed.
394 *
395 * @param b Parameters matrix that determines a linear system of equations.
396 * Provided matrix must have the same number of rows as provided input
397 * matrix for QR decomposition. Besides, each column on parameters matrix
398 * will represent a new system of equations, whose solution will be returned
399 * on appropriate column as an output of this method.
400 * @param roundingError Determines the amount of margin given to determine
401 * whether a matrix has full rank or not due to rounding errors. If not
402 * provided, by default rounding error is set to zero, but this value can be
403 * relaxed if needed.
404 * @param result Matrix containing solution of linear system of equations on
405 * each column for each column of provided parameters matrix b.
406 * @throws NotAvailableException Exception thrown if attempting to call this
407 * method before computing QR decomposition. To avoid this exception call
408 * decompose() method first.
409 * @throws WrongSizeException Exception thrown if attempting to call this
410 * method using an input matrix with less rows than columns; or if provided
411 * parameters matrix (b) does not have the same number of rows as input
412 * matrix being QR decomposed.
413 * @throws RankDeficientMatrixException Exception thrown if provided input
414 * matrix to be QR decomposed is rank deficient. In this case linear system
415 * of equations cannot be solved.
416 * @throws IllegalArgumentException Exception thrown if provided rounding
417 * error is lower than minimum allowed value (MIN_ROUND_ERROR).
418 * @see #decompose()
419 */
420 public void solve(final Matrix b, final double roundingError, final Matrix result)
421 throws NotAvailableException, WrongSizeException, RankDeficientMatrixException {
422
423 if (!isDecompositionAvailable()) {
424 throw new NotAvailableException();
425 }
426
427 final var rows = inputMatrix.getRows();
428 final var columns = inputMatrix.getColumns();
429 final var rowsB = b.getRows();
430 final var colsB = b.getColumns();
431 double sum;
432
433 if (rowsB != rows) {
434 throw new WrongSizeException();
435 }
436 if (roundingError < MIN_ROUND_ERROR) {
437 throw new IllegalArgumentException();
438 }
439 if (rows < columns) {
440 throw new WrongSizeException();
441 }
442 if (!isFullRank(roundingError)) {
443 throw new RankDeficientMatrixException();
444 }
445
446 // Compute Y = Q' * B
447 final var y = q.transposeAndReturnNew().multiplyAndReturnNew(b);
448
449 // resize result matrix if needed
450 if (result.getRows() != columns || result.getColumns() != colsB) {
451 result.resize(columns, colsB);
452 }
453
454 // Solve R * X = Y
455 // Each column of B will be a column of out (i.e. a solution of the
456 // linear system of equations)
457 for (var j2 = 0; j2 < colsB; j2++) {
458 // for overdetermined systems R has rows > columns, so we use only
459 // first columns rows, which contain upper diagonal data of R, the
460 // remaining rows of R are just zero.
461 for (var i = columns - 1; i >= 0; i--) {
462 sum = y.getElementAt(i, j2);
463 for (var j = i + 1; j < columns; j++) {
464 sum -= r.getElementAt(i, j) * result.getElementAt(j, j2);
465 }
466 result.setElementAt(i, j2, sum / r.getElementAt(i, i));
467 }
468 }
469 }
470
471 /**
472 * Solves a linear system of equations of the following form: A * X = B.
473 * Where A is the input matrix provided for QR decomposition, X is the
474 * solution to the system of equations, and B is the parameters vector/
475 * matrix.
476 * Note: This method can be reused for different b vectors/matrices without
477 * having to recompute QR decomposition on the same input matrix.
478 * Note: Provided b matrix must have the same number of rows as provided
479 * input matrix A, otherwise a WrongSizeException will be raised.
480 * Note: Provided input matrix "A" must have at least as many rows as columns,
481 * otherwise a WrongSizeException will be raised as well. For input matrices
482 * having a higher number of rows than columns, the system of equations will
483 * be overdetermined and the least squares solution will be found.
484 * Note: If provided input matrix A is rank deficient, a
485 * RankDeficientMatrixException will be thrown.
486 * Note: In order to execute this method, a QR decomposition must be
487 * available, otherwise a NotAvailableException will be raised. In order to
488 * avoid this exception call decompose() method first.
489 * Note: To solve the linear system of equations DEFAULT_ROUND_ERROR is used
490 *
491 * @param b Parameters matrix that determines a linear system of equations.
492 * Provided matrix must have the same number of rows as provided input
493 * matrix for QR decomposition. Besides, each column on parameters matrix
494 * will represent a new system of equations, whose solution will be returned
495 * on appropriate column as an output of this method.
496 * @return Matrix containing solution of linear system of equations on each
497 * column for each column of provided parameters matrix b.
498 * @throws NotAvailableException Exception thrown if attempting to call this
499 * method before computing QR decomposition. To avoid this exception call
500 * decompose() method first.
501 * @throws WrongSizeException Exception thrown if attempting to call this
502 * method using an input matrix with less rows than columns; or if provided
503 * parameters matrix (b) does not have the same number of rows as input
504 * matrix being QR decomposed.
505 * @throws RankDeficientMatrixException Exception thrown if provided input
506 * matrix to be QR decomposed is rank deficient. In this case linear system
507 * of equations cannot be solved.
508 * @see #decompose()
509 */
510 public Matrix solve(final Matrix b) throws NotAvailableException, WrongSizeException, RankDeficientMatrixException {
511 return solve(b, DEFAULT_ROUND_ERROR);
512 }
513
514 /**
515 * Solves a linear system of equations of the following form: A * X = B.
516 * Where A is the input matrix provided for QR decomposition, X is the
517 * solution to the system of equations, and B is the parameters vector/
518 * matrix.
519 * Note: This method can be reused for different b vectors/matrices without
520 * having to recompute QR decomposition on the same input matrix.
521 * Note: Provided b matrix must have the same number of rows as provided
522 * input matrix A, otherwise a WrongSizeException will be raised.
523 * Note: Provided input matrix "A" must have at least as many rows as columns,
524 * otherwise a WrongSizeException will be raised as well. For input matrices
525 * having a higher number of rows than columns, the system of equations will
526 * be overdetermined and the least squares solution will be found.
527 * Note: If provided input matrix A is rank deficient, a
528 * RankDeficientMatrixException will be thrown.
529 * Note: In order to execute this method, a QR decomposition must be
530 * available, otherwise a NotAvailableException will be raised. In order to
531 * avoid this exception call decompose() method first.
532 *
533 * @param b Parameters matrix that determines a linear system of equations.
534 * Provided matrix must have the same number of rows as provided input
535 * matrix for QR decomposition. Besides, each column on parameters matrix
536 * will represent a new system of equations, whose solution will be returned
537 * on appropriate column as an output of this method.
538 * @param roundingError Determines the amount of margin given to determine
539 * whether a matrix has full rank or not due to rounding errors. If not
540 * provided, by default rounding error is set to zero, but this value can be
541 * relaxed if needed.
542 * @return Matrix containing solution of linear system of equations on each
543 * column for each column of provided parameters matrix b.
544 * @throws NotAvailableException Exception thrown if attempting to call this
545 * method before computing QR decomposition. To avoid this exception call
546 * decompose() method first.
547 * @throws WrongSizeException Exception thrown if attempting to call this
548 * method using an input matrix with less rows than columns; or if provided
549 * parameters matrix (b) does not have the same number of rows as input
550 * matrix being QR decomposed.
551 * @throws RankDeficientMatrixException Exception thrown if provided input
552 * matrix to be QR decomposed is rank deficient. In this case linear system
553 * of equations cannot be solved.
554 * @throws IllegalArgumentException Exception thrown if provided rounding
555 * error is lower than minimum allowed value (MIN_ROUND_ERROR).
556 * @see #decompose()
557 */
558 public Matrix solve(final Matrix b, final double roundingError) throws NotAvailableException, WrongSizeException,
559 RankDeficientMatrixException {
560
561 if (!isDecompositionAvailable()) {
562 throw new NotAvailableException();
563 }
564
565 final var columns = inputMatrix.getColumns();
566 final var colsB = b.getColumns();
567 final var out = new Matrix(columns, colsB);
568 solve(b, roundingError, out);
569 return out;
570 }
571 }