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 java.util.ArrayList;
26 import java.util.Arrays;
27 import java.util.List;
28
29 import org.hipparchus.CalculusFieldElement;
30 import org.hipparchus.Field;
31 import org.hipparchus.FieldElement;
32 import org.hipparchus.exception.LocalizedCoreFormats;
33 import org.hipparchus.exception.MathIllegalArgumentException;
34 import org.hipparchus.exception.MathRuntimeException;
35 import org.hipparchus.exception.NullArgumentException;
36 import org.hipparchus.fraction.BigFraction;
37 import org.hipparchus.fraction.Fraction;
38 import org.hipparchus.util.FastMath;
39 import org.hipparchus.util.MathArrays;
40 import org.hipparchus.util.MathUtils;
41 import org.hipparchus.util.Precision;
42
43 /**
44 * A collection of static methods that operate on or return matrices.
45 *
46 */
47 public class MatrixUtils {
48
49 /**
50 * The default format for {@link RealMatrix} objects.
51 */
52 public static final RealMatrixFormat DEFAULT_FORMAT = RealMatrixFormat.getRealMatrixFormat();
53
54 /**
55 * A format for {@link RealMatrix} objects compatible with octave.
56 */
57 public static final RealMatrixFormat OCTAVE_FORMAT = new RealMatrixFormat("[", "]", "", "", "; ", ", ");
58
59 /** Pade coefficients required for the matrix exponential calculation. */
60 private static final double[] PADE_COEFFICIENTS_3 = {
61 120.0,
62 60.0,
63 12.0,
64 1.0
65 };
66
67 /** Pade coefficients required for the matrix exponential calculation. */
68 private static final double[] PADE_COEFFICIENTS_5 = {
69 30240.0,
70 15120.0,
71 3360.0,
72 420.0,
73 30.0,
74 1
75 };
76
77 /** Pade coefficients required for the matrix exponential calculation. */
78 private static final double[] PADE_COEFFICIENTS_7 = {
79 17297280.0,
80 8648640.0,
81 1995840.0,
82 277200.0,
83 25200.0,
84 1512.0,
85 56.0,
86 1.0
87 };
88
89 /** Pade coefficients required for the matrix exponential calculation. */
90 private static final double[] PADE_COEFFICIENTS_9 = {
91 17643225600.0,
92 8821612800.0,
93 2075673600.0,
94 302702400.0,
95 30270240.0,
96 2162160.0,
97 110880.0,
98 3960.0,
99 90.0,
100 1.0
101 };
102
103 /** Pade coefficients required for the matrix exponential calculation. */
104 private static final double[] PADE_COEFFICIENTS_13 = {
105 6.476475253248e+16,
106 3.238237626624e+16,
107 7.7717703038976e+15,
108 1.1873537964288e+15,
109 129060195264000.0,
110 10559470521600.0,
111 670442572800.0,
112 33522128640.0,
113 1323241920.0,
114 40840800.0,
115 960960.0,
116 16380.0,
117 182.0,
118 1.0
119 };
120
121 /**
122 * Private constructor.
123 */
124 private MatrixUtils() {
125 super();
126 }
127
128 /**
129 * Returns a {@link RealMatrix} with specified dimensions.
130 * <p>The type of matrix returned depends on the dimension. Below
131 * 2<sup>12</sup> elements (i.e. 4096 elements or 64×64 for a
132 * square matrix) which can be stored in a 32kB array, a {@link
133 * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
134 * BlockRealMatrix} instance is built.</p>
135 * <p>The matrix elements are all set to 0.0.</p>
136 * @param rows number of rows of the matrix
137 * @param columns number of columns of the matrix
138 * @return RealMatrix with specified dimensions
139 * @see #createRealMatrix(double[][])
140 */
141 public static RealMatrix createRealMatrix(final int rows, final int columns) {
142 return (rows * columns <= 4096) ?
143 new Array2DRowRealMatrix(rows, columns) : new BlockRealMatrix(rows, columns);
144 }
145
146 /**
147 * Returns a {@link FieldMatrix} with specified dimensions.
148 * <p>The type of matrix returned depends on the dimension. Below
149 * 2<sup>12</sup> elements (i.e. 4096 elements or 64×64 for a
150 * square matrix), a {@link FieldMatrix} instance is built. Above
151 * this threshold a {@link BlockFieldMatrix} instance is built.</p>
152 * <p>The matrix elements are all set to field.getZero().</p>
153 * @param <T> the type of the field elements
154 * @param field field to which the matrix elements belong
155 * @param rows number of rows of the matrix
156 * @param columns number of columns of the matrix
157 * @return FieldMatrix with specified dimensions
158 * @see #createFieldMatrix(FieldElement[][])
159 */
160 public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(final Field<T> field,
161 final int rows,
162 final int columns) {
163 return (rows * columns <= 4096) ?
164 new Array2DRowFieldMatrix<T>(field, rows, columns) : new BlockFieldMatrix<T>(field, rows, columns);
165 }
166
167 /**
168 * Returns a {@link RealMatrix} whose entries are the the values in the
169 * the input array.
170 * <p>The type of matrix returned depends on the dimension. Below
171 * 2<sup>12</sup> elements (i.e. 4096 elements or 64×64 for a
172 * square matrix) which can be stored in a 32kB array, a {@link
173 * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
174 * BlockRealMatrix} instance is built.</p>
175 * <p>The input array is copied, not referenced.</p>
176 *
177 * @param data input array
178 * @return RealMatrix containing the values of the array
179 * @throws org.hipparchus.exception.MathIllegalArgumentException
180 * if {@code data} is not rectangular (not all rows have the same length).
181 * @throws MathIllegalArgumentException if a row or column is empty.
182 * @throws NullArgumentException if either {@code data} or {@code data[0]}
183 * is {@code null}.
184 * @throws MathIllegalArgumentException if {@code data} is not rectangular.
185 * @see #createRealMatrix(int, int)
186 */
187 public static RealMatrix createRealMatrix(double[][] data)
188 throws MathIllegalArgumentException, NullArgumentException {
189 if (data == null ||
190 data[0] == null) {
191 throw new NullArgumentException();
192 }
193 return (data.length * data[0].length <= 4096) ?
194 new Array2DRowRealMatrix(data) : new BlockRealMatrix(data);
195 }
196
197 /**
198 * Returns a {@link FieldMatrix} whose entries are the the values in the
199 * the input array.
200 * <p>The type of matrix returned depends on the dimension. Below
201 * 2<sup>12</sup> elements (i.e. 4096 elements or 64×64 for a
202 * square matrix), a {@link FieldMatrix} instance is built. Above
203 * this threshold a {@link BlockFieldMatrix} instance is built.</p>
204 * <p>The input array is copied, not referenced.</p>
205 * @param <T> the type of the field elements
206 * @param data input array
207 * @return a matrix containing the values of the array.
208 * @throws org.hipparchus.exception.MathIllegalArgumentException
209 * if {@code data} is not rectangular (not all rows have the same length).
210 * @throws MathIllegalArgumentException if a row or column is empty.
211 * @throws NullArgumentException if either {@code data} or {@code data[0]}
212 * is {@code null}.
213 * @see #createFieldMatrix(Field, int, int)
214 */
215 public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(T[][] data)
216 throws MathIllegalArgumentException, NullArgumentException {
217 if (data == null ||
218 data[0] == null) {
219 throw new NullArgumentException();
220 }
221 return (data.length * data[0].length <= 4096) ?
222 new Array2DRowFieldMatrix<T>(data) : new BlockFieldMatrix<T>(data);
223 }
224
225 /**
226 * Returns <code>dimension x dimension</code> identity matrix.
227 *
228 * @param dimension dimension of identity matrix to generate
229 * @return identity matrix
230 * @throws IllegalArgumentException if dimension is not positive
231 */
232 public static RealMatrix createRealIdentityMatrix(int dimension) {
233 final RealMatrix m = createRealMatrix(dimension, dimension);
234 for (int i = 0; i < dimension; ++i) {
235 m.setEntry(i, i, 1.0);
236 }
237 return m;
238 }
239
240 /**
241 * Returns <code>dimension x dimension</code> identity matrix.
242 *
243 * @param <T> the type of the field elements
244 * @param field field to which the elements belong
245 * @param dimension dimension of identity matrix to generate
246 * @return identity matrix
247 * @throws IllegalArgumentException if dimension is not positive
248 */
249 public static <T extends FieldElement<T>> FieldMatrix<T>
250 createFieldIdentityMatrix(final Field<T> field, final int dimension) {
251 final T zero = field.getZero();
252 final T one = field.getOne();
253 final T[][] d = MathArrays.buildArray(field, dimension, dimension);
254 for (int row = 0; row < dimension; row++) {
255 final T[] dRow = d[row];
256 Arrays.fill(dRow, zero);
257 dRow[row] = one;
258 }
259 return new Array2DRowFieldMatrix<T>(field, d, false);
260 }
261
262 /**
263 * Returns a diagonal matrix with specified elements.
264 *
265 * @param diagonal diagonal elements of the matrix (the array elements
266 * will be copied)
267 * @return diagonal matrix
268 */
269 public static RealMatrix createRealDiagonalMatrix(final double[] diagonal) {
270 final RealMatrix m = createRealMatrix(diagonal.length, diagonal.length);
271 for (int i = 0; i < diagonal.length; ++i) {
272 m.setEntry(i, i, diagonal[i]);
273 }
274 return m;
275 }
276
277 /**
278 * Returns a diagonal matrix with specified elements.
279 *
280 * @param <T> the type of the field elements
281 * @param diagonal diagonal elements of the matrix (the array elements
282 * will be copied)
283 * @return diagonal matrix
284 */
285 public static <T extends FieldElement<T>> FieldMatrix<T>
286 createFieldDiagonalMatrix(final T[] diagonal) {
287 final FieldMatrix<T> m =
288 createFieldMatrix(diagonal[0].getField(), diagonal.length, diagonal.length);
289 for (int i = 0; i < diagonal.length; ++i) {
290 m.setEntry(i, i, diagonal[i]);
291 }
292 return m;
293 }
294
295 /**
296 * Creates a {@link RealVector} using the data from the input array.
297 *
298 * @param data the input data
299 * @return a data.length RealVector
300 * @throws MathIllegalArgumentException if {@code data} is empty.
301 * @throws NullArgumentException if {@code data} is {@code null}.
302 */
303 public static RealVector createRealVector(double[] data)
304 throws MathIllegalArgumentException, NullArgumentException {
305 if (data == null) {
306 throw new NullArgumentException();
307 }
308 return new ArrayRealVector(data, true);
309 }
310
311 /**
312 * Creates a {@link RealVector} with specified dimensions.
313 *
314 * @param dimension dimension of the vector
315 * @return a new vector
316 * @since 1.3
317 */
318 public static RealVector createRealVector(final int dimension) {
319 return new ArrayRealVector(new double[dimension]);
320 }
321
322 /**
323 * Creates a {@link FieldVector} using the data from the input array.
324 *
325 * @param <T> the type of the field elements
326 * @param data the input data
327 * @return a data.length FieldVector
328 * @throws MathIllegalArgumentException if {@code data} is empty.
329 * @throws NullArgumentException if {@code data} is {@code null}.
330 * @throws MathIllegalArgumentException if {@code data} has 0 elements
331 */
332 public static <T extends FieldElement<T>> FieldVector<T> createFieldVector(final T[] data)
333 throws MathIllegalArgumentException, NullArgumentException {
334 if (data == null) {
335 throw new NullArgumentException();
336 }
337 if (data.length == 0) {
338 throw new MathIllegalArgumentException(LocalizedCoreFormats.VECTOR_MUST_HAVE_AT_LEAST_ONE_ELEMENT);
339 }
340 return new ArrayFieldVector<T>(data[0].getField(), data, true);
341 }
342
343 /**
344 * Creates a {@link FieldVector} with specified dimensions.
345 *
346 * @param <T> the type of the field elements
347 * @param field field to which array elements belong
348 * @param dimension dimension of the vector
349 * @return a new vector
350 * @since 1.3
351 */
352 public static <T extends FieldElement<T>> FieldVector<T> createFieldVector(final Field<T> field, final int dimension) {
353 return new ArrayFieldVector<>(MathArrays.buildArray(field, dimension));
354 }
355
356 /**
357 * Create a row {@link RealMatrix} using the data from the input
358 * array.
359 *
360 * @param rowData the input row data
361 * @return a 1 x rowData.length RealMatrix
362 * @throws MathIllegalArgumentException if {@code rowData} is empty.
363 * @throws NullArgumentException if {@code rowData} is {@code null}.
364 */
365 public static RealMatrix createRowRealMatrix(double[] rowData)
366 throws MathIllegalArgumentException, NullArgumentException {
367 if (rowData == null) {
368 throw new NullArgumentException();
369 }
370 final int nCols = rowData.length;
371 final RealMatrix m = createRealMatrix(1, nCols);
372 for (int i = 0; i < nCols; ++i) {
373 m.setEntry(0, i, rowData[i]);
374 }
375 return m;
376 }
377
378 /**
379 * Create a row {@link FieldMatrix} using the data from the input
380 * array.
381 *
382 * @param <T> the type of the field elements
383 * @param rowData the input row data
384 * @return a 1 x rowData.length FieldMatrix
385 * @throws MathIllegalArgumentException if {@code rowData} is empty.
386 * @throws NullArgumentException if {@code rowData} is {@code null}.
387 */
388 public static <T extends FieldElement<T>> FieldMatrix<T>
389 createRowFieldMatrix(final T[] rowData)
390 throws MathIllegalArgumentException, NullArgumentException {
391 if (rowData == null) {
392 throw new NullArgumentException();
393 }
394 final int nCols = rowData.length;
395 if (nCols == 0) {
396 throw new MathIllegalArgumentException(LocalizedCoreFormats.AT_LEAST_ONE_COLUMN);
397 }
398 final FieldMatrix<T> m = createFieldMatrix(rowData[0].getField(), 1, nCols);
399 for (int i = 0; i < nCols; ++i) {
400 m.setEntry(0, i, rowData[i]);
401 }
402 return m;
403 }
404
405 /**
406 * Creates a column {@link RealMatrix} using the data from the input
407 * array.
408 *
409 * @param columnData the input column data
410 * @return a columnData x 1 RealMatrix
411 * @throws MathIllegalArgumentException if {@code columnData} is empty.
412 * @throws NullArgumentException if {@code columnData} is {@code null}.
413 */
414 public static RealMatrix createColumnRealMatrix(double[] columnData)
415 throws MathIllegalArgumentException, NullArgumentException {
416 if (columnData == null) {
417 throw new NullArgumentException();
418 }
419 final int nRows = columnData.length;
420 final RealMatrix m = createRealMatrix(nRows, 1);
421 for (int i = 0; i < nRows; ++i) {
422 m.setEntry(i, 0, columnData[i]);
423 }
424 return m;
425 }
426
427 /**
428 * Creates a column {@link FieldMatrix} using the data from the input
429 * array.
430 *
431 * @param <T> the type of the field elements
432 * @param columnData the input column data
433 * @return a columnData x 1 FieldMatrix
434 * @throws MathIllegalArgumentException if {@code data} is empty.
435 * @throws NullArgumentException if {@code columnData} is {@code null}.
436 */
437 public static <T extends FieldElement<T>> FieldMatrix<T>
438 createColumnFieldMatrix(final T[] columnData)
439 throws MathIllegalArgumentException, NullArgumentException {
440 if (columnData == null) {
441 throw new NullArgumentException();
442 }
443 final int nRows = columnData.length;
444 if (nRows == 0) {
445 throw new MathIllegalArgumentException(LocalizedCoreFormats.AT_LEAST_ONE_ROW);
446 }
447 final FieldMatrix<T> m = createFieldMatrix(columnData[0].getField(), nRows, 1);
448 for (int i = 0; i < nRows; ++i) {
449 m.setEntry(i, 0, columnData[i]);
450 }
451 return m;
452 }
453
454 /**
455 * Checks whether a matrix is symmetric, within a given relative tolerance.
456 *
457 * @param matrix Matrix to check.
458 * @param relativeTolerance Tolerance of the symmetry check.
459 * @param raiseException If {@code true}, an exception will be raised if
460 * the matrix is not symmetric.
461 * @return {@code true} if {@code matrix} is symmetric.
462 * @throws MathIllegalArgumentException if the matrix is not square.
463 * @throws MathIllegalArgumentException if the matrix is not symmetric.
464 */
465 private static boolean isSymmetricInternal(RealMatrix matrix,
466 double relativeTolerance,
467 boolean raiseException) {
468 final int rows = matrix.getRowDimension();
469 if (rows != matrix.getColumnDimension()) {
470 if (raiseException) {
471 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
472 rows, matrix.getColumnDimension());
473 } else {
474 return false;
475 }
476 }
477 for (int i = 0; i < rows; i++) {
478 for (int j = i + 1; j < rows; j++) {
479 final double mij = matrix.getEntry(i, j);
480 final double mji = matrix.getEntry(j, i);
481 if (FastMath.abs(mij - mji) >
482 FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * relativeTolerance) {
483 if (raiseException) {
484 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SYMMETRIC_MATRIX,
485 i, j, relativeTolerance);
486 } else {
487 return false;
488 }
489 }
490 }
491 }
492 return true;
493 }
494
495 /**
496 * Checks whether a matrix is symmetric.
497 *
498 * @param matrix Matrix to check.
499 * @param eps Relative tolerance.
500 * @throws MathIllegalArgumentException if the matrix is not square.
501 * @throws MathIllegalArgumentException if the matrix is not symmetric.
502 */
503 public static void checkSymmetric(RealMatrix matrix,
504 double eps) {
505 isSymmetricInternal(matrix, eps, true);
506 }
507
508 /**
509 * Checks whether a matrix is symmetric.
510 *
511 * @param matrix Matrix to check.
512 * @param eps Relative tolerance.
513 * @return {@code true} if {@code matrix} is symmetric.
514 */
515 public static boolean isSymmetric(RealMatrix matrix,
516 double eps) {
517 return isSymmetricInternal(matrix, eps, false);
518 }
519
520 /**
521 * Check if matrix indices are valid.
522 *
523 * @param m Matrix.
524 * @param row Row index to check.
525 * @param column Column index to check.
526 * @throws MathIllegalArgumentException if {@code row} or {@code column} is not
527 * a valid index.
528 */
529 public static void checkMatrixIndex(final AnyMatrix m,
530 final int row, final int column)
531 throws MathIllegalArgumentException {
532 checkRowIndex(m, row);
533 checkColumnIndex(m, column);
534 }
535
536 /**
537 * Check if a row index is valid.
538 *
539 * @param m Matrix.
540 * @param row Row index to check.
541 * @throws MathIllegalArgumentException if {@code row} is not a valid index.
542 */
543 public static void checkRowIndex(final AnyMatrix m, final int row)
544 throws MathIllegalArgumentException {
545 if (row < 0 ||
546 row >= m.getRowDimension()) {
547 throw new MathIllegalArgumentException(LocalizedCoreFormats.ROW_INDEX,
548 row, 0, m.getRowDimension() - 1);
549 }
550 }
551
552 /**
553 * Check if a column index is valid.
554 *
555 * @param m Matrix.
556 * @param column Column index to check.
557 * @throws MathIllegalArgumentException if {@code column} is not a valid index.
558 */
559 public static void checkColumnIndex(final AnyMatrix m, final int column)
560 throws MathIllegalArgumentException {
561 if (column < 0 || column >= m.getColumnDimension()) {
562 throw new MathIllegalArgumentException(LocalizedCoreFormats.COLUMN_INDEX,
563 column, 0, m.getColumnDimension() - 1);
564 }
565 }
566
567 /**
568 * Check if submatrix ranges indices are valid.
569 * Rows and columns are indicated counting from 0 to {@code n - 1}.
570 *
571 * @param m Matrix.
572 * @param startRow Initial row index.
573 * @param endRow Final row index.
574 * @param startColumn Initial column index.
575 * @param endColumn Final column index.
576 * @throws MathIllegalArgumentException if the indices are invalid.
577 * @throws MathIllegalArgumentException if {@code endRow < startRow} or
578 * {@code endColumn < startColumn}.
579 */
580 public static void checkSubMatrixIndex(final AnyMatrix m,
581 final int startRow, final int endRow,
582 final int startColumn, final int endColumn)
583 throws MathIllegalArgumentException {
584 checkRowIndex(m, startRow);
585 checkRowIndex(m, endRow);
586 if (endRow < startRow) {
587 throw new MathIllegalArgumentException(LocalizedCoreFormats.INITIAL_ROW_AFTER_FINAL_ROW,
588 endRow, startRow, false);
589 }
590
591 checkColumnIndex(m, startColumn);
592 checkColumnIndex(m, endColumn);
593 if (endColumn < startColumn) {
594 throw new MathIllegalArgumentException(LocalizedCoreFormats.INITIAL_COLUMN_AFTER_FINAL_COLUMN,
595 endColumn, startColumn, false);
596 }
597
598
599 }
600
601 /**
602 * Check if submatrix ranges indices are valid.
603 * Rows and columns are indicated counting from 0 to n-1.
604 *
605 * @param m Matrix.
606 * @param selectedRows Array of row indices.
607 * @param selectedColumns Array of column indices.
608 * @throws NullArgumentException if {@code selectedRows} or
609 * {@code selectedColumns} are {@code null}.
610 * @throws MathIllegalArgumentException if the row or column selections are empty (zero
611 * length).
612 * @throws MathIllegalArgumentException if row or column selections are not valid.
613 */
614 public static void checkSubMatrixIndex(final AnyMatrix m,
615 final int[] selectedRows,
616 final int[] selectedColumns)
617 throws MathIllegalArgumentException, NullArgumentException {
618 if (selectedRows == null) {
619 throw new NullArgumentException();
620 }
621 if (selectedColumns == null) {
622 throw new NullArgumentException();
623 }
624 if (selectedRows.length == 0) {
625 throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_SELECTED_ROW_INDEX_ARRAY);
626 }
627 if (selectedColumns.length == 0) {
628 throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_SELECTED_COLUMN_INDEX_ARRAY);
629 }
630
631 for (final int row : selectedRows) {
632 checkRowIndex(m, row);
633 }
634 for (final int column : selectedColumns) {
635 checkColumnIndex(m, column);
636 }
637 }
638
639 /**
640 * Check if matrices are addition compatible.
641 *
642 * @param left Left hand side matrix.
643 * @param right Right hand side matrix.
644 * @throws MathIllegalArgumentException if the matrices are not addition
645 * compatible.
646 */
647 public static void checkAdditionCompatible(final AnyMatrix left, final AnyMatrix right)
648 throws MathIllegalArgumentException {
649 if ((left.getRowDimension() != right.getRowDimension()) ||
650 (left.getColumnDimension() != right.getColumnDimension())) {
651 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
652 left.getRowDimension(), left.getColumnDimension(),
653 right.getRowDimension(), right.getColumnDimension());
654 }
655 }
656
657 /**
658 * Check if matrices are subtraction compatible
659 *
660 * @param left Left hand side matrix.
661 * @param right Right hand side matrix.
662 * @throws MathIllegalArgumentException if the matrices are not addition
663 * compatible.
664 */
665 public static void checkSubtractionCompatible(final AnyMatrix left, final AnyMatrix right)
666 throws MathIllegalArgumentException {
667 if ((left.getRowDimension() != right.getRowDimension()) ||
668 (left.getColumnDimension() != right.getColumnDimension())) {
669 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
670 left.getRowDimension(), left.getColumnDimension(),
671 right.getRowDimension(), right.getColumnDimension());
672 }
673 }
674
675 /**
676 * Check if matrices are multiplication compatible
677 *
678 * @param left Left hand side matrix.
679 * @param right Right hand side matrix.
680 * @throws MathIllegalArgumentException if matrices are not multiplication
681 * compatible.
682 */
683 public static void checkMultiplicationCompatible(final AnyMatrix left, final AnyMatrix right)
684 throws MathIllegalArgumentException {
685 if (left.getColumnDimension() != right.getRowDimension()) {
686 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
687 left.getColumnDimension(), right.getRowDimension());
688 }
689 }
690
691 /**
692 * Check if matrices have the same number of columns.
693 *
694 * @param left Left hand side matrix.
695 * @param right Right hand side matrix.
696 * @throws MathIllegalArgumentException if matrices don't have the same number of columns.
697 * @since 1.3
698 */
699 public static void checkSameColumnDimension(final AnyMatrix left, final AnyMatrix right)
700 throws MathIllegalArgumentException {
701 if (left.getColumnDimension() != right.getColumnDimension()) {
702 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
703 left.getColumnDimension(), right.getColumnDimension());
704 }
705 }
706
707 /**
708 * Check if matrices have the same number of rows.
709 *
710 * @param left Left hand side matrix.
711 * @param right Right hand side matrix.
712 * @throws MathIllegalArgumentException if matrices don't have the same number of rows.
713 * @since 1.3
714 */
715 public static void checkSameRowDimension(final AnyMatrix left, final AnyMatrix right)
716 throws MathIllegalArgumentException {
717 if (left.getRowDimension() != right.getRowDimension()) {
718 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
719 left.getRowDimension(), right.getRowDimension());
720 }
721 }
722
723 /**
724 * Convert a {@link FieldMatrix}/{@link Fraction} matrix to a {@link RealMatrix}.
725 * @param m Matrix to convert.
726 * @return the converted matrix.
727 */
728 public static Array2DRowRealMatrix fractionMatrixToRealMatrix(final FieldMatrix<Fraction> m) {
729 final FractionMatrixConverter converter = new FractionMatrixConverter();
730 m.walkInOptimizedOrder(converter);
731 return converter.getConvertedMatrix();
732 }
733
734 /** Converter for {@link FieldMatrix}/{@link Fraction}. */
735 private static class FractionMatrixConverter extends DefaultFieldMatrixPreservingVisitor<Fraction> {
736 /** Converted array. */
737 private double[][] data;
738 /** Simple constructor. */
739 FractionMatrixConverter() {
740 super(Fraction.ZERO);
741 }
742
743 /** {@inheritDoc} */
744 @Override
745 public void start(int rows, int columns,
746 int startRow, int endRow, int startColumn, int endColumn) {
747 data = new double[rows][columns];
748 }
749
750 /** {@inheritDoc} */
751 @Override
752 public void visit(int row, int column, Fraction value) {
753 data[row][column] = value.doubleValue();
754 }
755
756 /**
757 * Get the converted matrix.
758 *
759 * @return the converted matrix.
760 */
761 Array2DRowRealMatrix getConvertedMatrix() {
762 return new Array2DRowRealMatrix(data, false);
763 }
764
765 }
766
767 /**
768 * Convert a {@link FieldMatrix}/{@link BigFraction} matrix to a {@link RealMatrix}.
769 *
770 * @param m Matrix to convert.
771 * @return the converted matrix.
772 */
773 public static Array2DRowRealMatrix bigFractionMatrixToRealMatrix(final FieldMatrix<BigFraction> m) {
774 final BigFractionMatrixConverter converter = new BigFractionMatrixConverter();
775 m.walkInOptimizedOrder(converter);
776 return converter.getConvertedMatrix();
777 }
778
779 /** Converter for {@link FieldMatrix}/{@link BigFraction}. */
780 private static class BigFractionMatrixConverter extends DefaultFieldMatrixPreservingVisitor<BigFraction> {
781 /** Converted array. */
782 private double[][] data;
783 /** Simple constructor. */
784 BigFractionMatrixConverter() {
785 super(BigFraction.ZERO);
786 }
787
788 /** {@inheritDoc} */
789 @Override
790 public void start(int rows, int columns,
791 int startRow, int endRow, int startColumn, int endColumn) {
792 data = new double[rows][columns];
793 }
794
795 /** {@inheritDoc} */
796 @Override
797 public void visit(int row, int column, BigFraction value) {
798 data[row][column] = value.doubleValue();
799 }
800
801 /**
802 * Get the converted matrix.
803 *
804 * @return the converted matrix.
805 */
806 Array2DRowRealMatrix getConvertedMatrix() {
807 return new Array2DRowRealMatrix(data, false);
808 }
809 }
810
811 /**Solve a system of composed of a Lower Triangular Matrix
812 * {@link RealMatrix}.
813 * <p>
814 * This method is called to solve systems of equations which are
815 * of the lower triangular form. The matrix {@link RealMatrix}
816 * is assumed, though not checked, to be in lower triangular form.
817 * The vector {@link RealVector} is overwritten with the solution.
818 * The matrix is checked that it is square and its dimensions match
819 * the length of the vector.
820 * </p>
821 * @param rm RealMatrix which is lower triangular
822 * @param b RealVector this is overwritten
823 * @throws MathIllegalArgumentException if the matrix and vector are not
824 * conformable
825 * @throws MathIllegalArgumentException if the matrix {@code rm} is not square
826 * @throws MathRuntimeException if the absolute value of one of the diagonal
827 * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
828 */
829 public static void solveLowerTriangularSystem(RealMatrix rm, RealVector b)
830 throws MathIllegalArgumentException, MathRuntimeException {
831 if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
832 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
833 (rm == null) ? 0 : rm.getRowDimension(),
834 (b == null) ? 0 : b.getDimension());
835 }
836 if( rm.getColumnDimension() != rm.getRowDimension() ){
837 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
838 rm.getRowDimension(), rm.getColumnDimension());
839 }
840 int rows = rm.getRowDimension();
841 for( int i = 0 ; i < rows ; i++ ){
842 double diag = rm.getEntry(i, i);
843 if( FastMath.abs(diag) < Precision.SAFE_MIN ){
844 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
845 }
846 double bi = b.getEntry(i)/diag;
847 b.setEntry(i, bi );
848 for( int j = i+1; j< rows; j++ ){
849 b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i) );
850 }
851 }
852 }
853
854 /** Solver a system composed of an Upper Triangular Matrix
855 * {@link RealMatrix}.
856 * <p>
857 * This method is called to solve systems of equations which are
858 * of the lower triangular form. The matrix {@link RealMatrix}
859 * is assumed, though not checked, to be in upper triangular form.
860 * The vector {@link RealVector} is overwritten with the solution.
861 * The matrix is checked that it is square and its dimensions match
862 * the length of the vector.
863 * </p>
864 * @param rm RealMatrix which is upper triangular
865 * @param b RealVector this is overwritten
866 * @throws MathIllegalArgumentException if the matrix and vector are not
867 * conformable
868 * @throws MathIllegalArgumentException if the matrix {@code rm} is not
869 * square
870 * @throws MathRuntimeException if the absolute value of one of the diagonal
871 * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
872 */
873 public static void solveUpperTriangularSystem(RealMatrix rm, RealVector b)
874 throws MathIllegalArgumentException, MathRuntimeException {
875 if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
876 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
877 (rm == null) ? 0 : rm.getRowDimension(),
878 (b == null) ? 0 : b.getDimension());
879 }
880 if( rm.getColumnDimension() != rm.getRowDimension() ){
881 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
882 rm.getRowDimension(), rm.getColumnDimension());
883 }
884 int rows = rm.getRowDimension();
885 for( int i = rows-1 ; i >-1 ; i-- ){
886 double diag = rm.getEntry(i, i);
887 if( FastMath.abs(diag) < Precision.SAFE_MIN ){
888 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
889 }
890 double bi = b.getEntry(i)/diag;
891 b.setEntry(i, bi );
892 for( int j = i-1; j>-1; j-- ){
893 b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i) );
894 }
895 }
896 }
897
898 /**
899 * Computes the inverse of the given matrix by splitting it into
900 * 4 sub-matrices.
901 *
902 * @param m Matrix whose inverse must be computed.
903 * @param splitIndex Index that determines the "split" line and
904 * column.
905 * The element corresponding to this index will part of the
906 * upper-left sub-matrix.
907 * @return the inverse of {@code m}.
908 * @throws MathIllegalArgumentException if {@code m} is not square.
909 */
910 public static RealMatrix blockInverse(RealMatrix m,
911 int splitIndex) {
912 final int n = m.getRowDimension();
913 if (m.getColumnDimension() != n) {
914 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
915 m.getRowDimension(), m.getColumnDimension());
916 }
917
918 final int splitIndex1 = splitIndex + 1;
919
920 final RealMatrix a = m.getSubMatrix(0, splitIndex, 0, splitIndex);
921 final RealMatrix b = m.getSubMatrix(0, splitIndex, splitIndex1, n - 1);
922 final RealMatrix c = m.getSubMatrix(splitIndex1, n - 1, 0, splitIndex);
923 final RealMatrix d = m.getSubMatrix(splitIndex1, n - 1, splitIndex1, n - 1);
924
925 final SingularValueDecomposition aDec = new SingularValueDecomposition(a);
926 final DecompositionSolver aSolver = aDec.getSolver();
927 if (!aSolver.isNonSingular()) {
928 throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
929 }
930 final RealMatrix aInv = aSolver.getInverse();
931
932 final SingularValueDecomposition dDec = new SingularValueDecomposition(d);
933 final DecompositionSolver dSolver = dDec.getSolver();
934 if (!dSolver.isNonSingular()) {
935 throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
936 }
937 final RealMatrix dInv = dSolver.getInverse();
938
939 final RealMatrix tmp1 = a.subtract(b.multiply(dInv).multiply(c));
940 final SingularValueDecomposition tmp1Dec = new SingularValueDecomposition(tmp1);
941 final DecompositionSolver tmp1Solver = tmp1Dec.getSolver();
942 if (!tmp1Solver.isNonSingular()) {
943 throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
944 }
945 final RealMatrix result00 = tmp1Solver.getInverse();
946
947 final RealMatrix tmp2 = d.subtract(c.multiply(aInv).multiply(b));
948 final SingularValueDecomposition tmp2Dec = new SingularValueDecomposition(tmp2);
949 final DecompositionSolver tmp2Solver = tmp2Dec.getSolver();
950 if (!tmp2Solver.isNonSingular()) {
951 throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
952 }
953 final RealMatrix result11 = tmp2Solver.getInverse();
954
955 final RealMatrix result01 = aInv.multiply(b).multiply(result11).scalarMultiply(-1);
956 final RealMatrix result10 = dInv.multiply(c).multiply(result00).scalarMultiply(-1);
957
958 final RealMatrix result = new Array2DRowRealMatrix(n, n);
959 result.setSubMatrix(result00.getData(), 0, 0);
960 result.setSubMatrix(result01.getData(), 0, splitIndex1);
961 result.setSubMatrix(result10.getData(), splitIndex1, 0);
962 result.setSubMatrix(result11.getData(), splitIndex1, splitIndex1);
963
964 return result;
965 }
966
967 /**
968 * Computes the inverse of the given matrix.
969 * <p>
970 * By default, the inverse of the matrix is computed using the QR-decomposition,
971 * unless a more efficient method can be determined for the input matrix.
972 * <p>
973 * Note: this method will use a singularity threshold of 0,
974 * use {@link #inverse(RealMatrix, double)} if a different threshold is needed.
975 *
976 * @param matrix Matrix whose inverse shall be computed
977 * @return the inverse of {@code matrix}
978 * @throws NullArgumentException if {@code matrix} is {@code null}
979 * @throws MathIllegalArgumentException if m is singular
980 * @throws MathIllegalArgumentException if matrix is not square
981 */
982 public static RealMatrix inverse(RealMatrix matrix)
983 throws MathIllegalArgumentException, NullArgumentException {
984 return inverse(matrix, 0);
985 }
986
987 /**
988 * Computes the inverse of the given matrix.
989 * <p>
990 * By default, the inverse of the matrix is computed using the QR-decomposition,
991 * unless a more efficient method can be determined for the input matrix.
992 *
993 * @param matrix Matrix whose inverse shall be computed
994 * @param threshold Singularity threshold
995 * @return the inverse of {@code m}
996 * @throws NullArgumentException if {@code matrix} is {@code null}
997 * @throws MathIllegalArgumentException if matrix is singular
998 * @throws MathIllegalArgumentException if matrix is not square
999 */
1000 public static RealMatrix inverse(RealMatrix matrix, double threshold)
1001 throws MathIllegalArgumentException, NullArgumentException {
1002
1003 MathUtils.checkNotNull(matrix);
1004
1005 if (!matrix.isSquare()) {
1006 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
1007 matrix.getRowDimension(), matrix.getColumnDimension());
1008 }
1009
1010 if (matrix instanceof DiagonalMatrix) {
1011 return ((DiagonalMatrix) matrix).inverse(threshold);
1012 } else {
1013 QRDecomposition decomposition = new QRDecomposition(matrix, threshold);
1014 return decomposition.getSolver().getInverse();
1015 }
1016 }
1017
1018 /**
1019 * Computes the <a href="https://mathworld.wolfram.com/MatrixExponential.html">
1020 * matrix exponential</a> of the given matrix.
1021 *
1022 * The algorithm implementation follows the Pade approximant method of
1023 * <p>Higham, Nicholas J. “The Scaling and Squaring Method for the Matrix Exponential
1024 * Revisited.” SIAM Journal on Matrix Analysis and Applications 26, no. 4 (January 2005): 1179–93.</p>
1025 *
1026 * @param rm RealMatrix whose inverse shall be computed
1027 * @return The inverse of {@code rm}
1028 * @throws MathIllegalArgumentException if matrix is not square
1029 */
1030 public static RealMatrix matrixExponential(final RealMatrix rm) {
1031
1032 // Check that the input matrix is square
1033 if (!rm.isSquare()) {
1034 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
1035 rm.getRowDimension(), rm.getColumnDimension());
1036 }
1037
1038 // Copy input matrix
1039 int dim = rm.getRowDimension();
1040 final RealMatrix identity = MatrixUtils.createRealIdentityMatrix(dim);
1041 RealMatrix scaledMatrix = rm.copy();
1042
1043 // Select pade degree required
1044 final double l1Norm = rm.getNorm1();
1045 double[] padeCoefficients;
1046 int squaringCount = 0;
1047
1048 if (l1Norm < 1.495585217958292e-2) {
1049 padeCoefficients = PADE_COEFFICIENTS_3;
1050 } else if (l1Norm < 2.539398330063230e-1) {
1051 padeCoefficients = PADE_COEFFICIENTS_5;
1052 } else if (l1Norm < 9.504178996162932e-1) {
1053 padeCoefficients = PADE_COEFFICIENTS_7;
1054 } else if (l1Norm < 2.097847961257068) {
1055 padeCoefficients = PADE_COEFFICIENTS_9;
1056 } else {
1057 padeCoefficients = PADE_COEFFICIENTS_13;
1058
1059 // Calculate scaling factor
1060 final double normScale = 5.371920351148152;
1061 squaringCount = Math.max(0, Math.getExponent(l1Norm / normScale));
1062
1063 // Scale matrix by power of 2
1064 final int finalSquaringCount = squaringCount;
1065 scaledMatrix.mapToSelf(x -> Math.scalb(x, -finalSquaringCount));
1066 }
1067
1068 // Calculate U and V using Horner
1069 // See Golub, Gene H., and Charles F. van Loan. Matrix Computations. 4th ed.
1070 // John Hopkins University Press, 2013. pages 530/531
1071 final RealMatrix scaledMatrix2 = scaledMatrix.multiply(scaledMatrix);
1072 final int coeffLength = padeCoefficients.length;
1073
1074 // Calculate V
1075 RealMatrix padeV = MatrixUtils.createRealMatrix(dim, dim);
1076 for (int i = coeffLength - 1; i > 1; i -= 2) {
1077 padeV = scaledMatrix2.multiply(padeV.add(identity.scalarMultiply(padeCoefficients[i])));
1078 }
1079 padeV = scaledMatrix.multiply(padeV.add(identity.scalarMultiply(padeCoefficients[1])));
1080
1081 // Calculate U
1082 RealMatrix padeU = MatrixUtils.createRealMatrix(dim, dim);
1083 for (int i = coeffLength - 2; i > 1; i -= 2) {
1084 padeU = scaledMatrix2.multiply(padeU.add(identity.scalarMultiply(padeCoefficients[i])));
1085 }
1086 padeU = padeU.add(identity.scalarMultiply(padeCoefficients[0]));
1087
1088 // Calculate pade approximate by solving (U-V) F = (U+V) for F
1089 RealMatrix padeNumer = padeU.add(padeV);
1090 RealMatrix padeDenom = padeU.subtract(padeV);
1091
1092 // Calculate the matrix ratio
1093 QRDecomposition decomposition = new QRDecomposition(padeDenom);
1094 RealMatrix result = decomposition.getSolver().solve(padeNumer);
1095
1096 // Repeated squaring if matrix was scaled
1097 for (int i = 0; i < squaringCount; i++) {
1098 result = result.multiply(result);
1099 }
1100
1101 return result;
1102 }
1103
1104 /** Orthonormalize a list of vectors.
1105 * <p>
1106 * Orthonormalization is performed by using the Modified Gram-Schmidt process.
1107 * </p>
1108 * @param independent list of independent vectors
1109 * @param threshold projected vectors with a norm less than or equal to this threshold
1110 * are considered to have zero norm, hence the vectors they come from are not independent from
1111 * previous vectors
1112 * @param handler handler for dependent vectors
1113 * @return orthonormal basis having the same span as {@code independent}
1114 * @since 2.1
1115 */
1116 public static List<RealVector> orthonormalize(final List<RealVector> independent,
1117 final double threshold, final DependentVectorsHandler handler) {
1118
1119 // create separate list
1120 final List<RealVector> basis = new ArrayList<>(independent);
1121
1122 // loop over basis vectors
1123 int index = 0;
1124 while (index < basis.size()) {
1125
1126 // check dependency
1127 final RealVector vi = basis.get(index);
1128 final double norm = vi.getNorm();
1129 if (norm <= threshold) {
1130 // the current vector is dependent from the previous ones
1131 index = handler.manageDependent(index, basis);
1132 } else {
1133
1134 // normalize basis vector in place
1135 vi.mapDivideToSelf(vi.getNorm());
1136
1137 // project remaining vectors in place
1138 for (int j = index + 1; j < basis.size(); ++j) {
1139 final RealVector vj = basis.get(j);
1140 final double dot = vi.dotProduct(vj);
1141 for (int k = 0; k < vj.getDimension(); ++k) {
1142 vj.setEntry(k, vj.getEntry(k) - dot * vi.getEntry(k));
1143 }
1144 }
1145
1146 ++index;
1147
1148 }
1149
1150 }
1151
1152 return basis;
1153
1154 }
1155
1156 /** Orthonormalize a list of vectors.
1157 * <p>
1158 * Orthonormalization is performed by using the Modified Gram-Schmidt process.
1159 * </p>
1160 * @param <T> type of the field elements
1161 * @param independent list of independent vectors
1162 * @param threshold projected vectors with a norm less than or equal to this threshold
1163 * are considered to have zero norm, hence the vectors they come from are not independent from
1164 * previous vectors
1165 * @param field type of the files elements
1166 * @param handler handler for dependent vectors
1167 * @return orthonormal basis having the same span as {@code independent}
1168 * @since 2.1
1169 */
1170 public static <T extends CalculusFieldElement<T>> List<FieldVector<T>> orthonormalize(final Field<T> field,
1171 final List<FieldVector<T>> independent,
1172 final T threshold,
1173 final DependentVectorsHandler handler) {
1174
1175 // create separate list
1176 final List<FieldVector<T>> basis = new ArrayList<>(independent);
1177
1178 // loop over basis vectors
1179 int index = 0;
1180 while (index < basis.size()) {
1181
1182 // check dependency
1183 final FieldVector<T> vi = basis.get(index);
1184 final T norm = vi.dotProduct(vi).sqrt();
1185 if (norm.subtract(threshold).getReal() <= 0) {
1186 // the current vector is dependent from the previous ones
1187 index = handler.manageDependent(field, index, basis);
1188 } else {
1189
1190 // normalize basis vector in place
1191 vi.mapDivideToSelf(norm);
1192
1193 // project remaining vectors in place
1194 for (int j = index + 1; j < basis.size(); ++j) {
1195 final FieldVector<T> vj = basis.get(j);
1196 final T dot = vi.dotProduct(vj);
1197 for (int k = 0; k < vj.getDimension(); ++k) {
1198 vj.setEntry(k, vj.getEntry(k).subtract(dot.multiply(vi.getEntry(k))));
1199 }
1200 }
1201
1202 ++index;
1203
1204 }
1205 }
1206
1207 return basis;
1208
1209 }
1210
1211 }