Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
algebra.h
Go to the documentation of this file.
1
18
19#ifndef THEORETICA_ALGEBRA_H
20#define THEORETICA_ALGEBRA_H
21
22#include "../complex/complex_types.h"
23#include "../core/core_traits.h"
24#include "../core/error.h"
25
26
27namespace theoretica {
28
29
31 namespace algebra {
32
33
34 // Operations involving one matrix or vector
35
36
43 template<typename Matrix>
44 inline Matrix& mat_error(Matrix& m) {
45
46 using Type = matrix_element_t<Matrix>;
47
48 for (unsigned int i = 0; i < m.rows(); ++i)
49 for (unsigned int j = 0; j < m.cols(); ++j)
50 m(i, j) = (Type) (i == j ? nan() : 0);
51
52 return m;
53 }
54
55
61 template<typename Vector>
62 inline Vector& vec_error(Vector& v) {
63
64 using Type = vector_element_t<Vector>;
65
66 for (unsigned int i = 0; i < v.size(); ++i)
67 v[i] = Type(nan());
68
69 return v;
70 }
71
72
76 template<typename Matrix>
77 inline Matrix& make_identity(Matrix& m) {
78
79 using Type = matrix_element_t<Matrix>;
80
81 for (unsigned int i = 0; i < m.rows(); ++i)
82 for (unsigned int j = 0; j < m.cols(); ++j)
83 m(i, j) = (Type) (i == j ? 1 : 0);
84
85 return m;
86 }
87
88
92 template<typename Matrix>
93 inline Matrix& mat_zeroes(Matrix& m) {
94
95 using Type = matrix_element_t<Matrix>;
96
97 for (unsigned int i = 0; i < m.rows(); ++i)
98 for (unsigned int j = 0; j < m.cols(); ++j)
99 m(i, j) = (Type) 0;
100
101 return m;
102 }
103
104
108 template<typename Vector>
109 inline Vector& vec_zeroes(Vector& v) {
110
111 using Type = vector_element_t<Vector>;
112
113 for (unsigned int i = 0; i < v.size(); ++i)
114 v[i] = Type(0.0);
115
116 return v;
117 }
118
119
124 template<typename Matrix1, typename Matrix2>
125 inline Matrix1& mat_copy(Matrix1& dest, const Matrix2& src) {
126
127 dest.resize(src.rows(), src.cols());
128
129 for (unsigned int i = 0; i < src.rows(); ++i)
130 for (unsigned int j = 0; j < src.cols(); ++j)
131 dest(i, j) = src(i, j);
132
133 return dest;
134 }
135
136
142 template<typename Vector1, typename Vector2>
143 inline Vector1& vec_copy(Vector1& dest, const Vector2& src) {
144
145 dest.resize(src.size());
146
147 for (unsigned int i = 0; i < src.size(); ++i)
148 dest[i] = src[i];
149
150 return dest;
151 }
152
153
160 template<typename Matrix>
161 inline Matrix& mat_swap_rows(Matrix& A, unsigned int row1, unsigned int row2) {
162
163 using Type = matrix_element_t<Matrix>;
164
165 if (row1 >= A.rows()) {
166 TH_MATH_ERROR("algebra::mat_swap_rows", row1, INVALID_ARGUMENT);
167 return mat_error(A);
168 }
169
170 if (row2 >= A.rows()) {
171 TH_MATH_ERROR("algebra::mat_swap_rows", row2, INVALID_ARGUMENT);
172 return mat_error(A);
173 }
174
175 if (row1 == row2)
176 return A;
177
178 for (unsigned int j = 0; j < A.cols(); ++j) {
179
180 const Type tmp = A(row1, j);
181
182 A(row1, j) = A(row2, j);
183 A(row2, j) = tmp;
184 }
185
186 return A;
187 }
188
189
196 template<typename Matrix>
197 inline Matrix& mat_swap_cols(Matrix& A, unsigned int col1, unsigned int col2) {
198
199 using Type = matrix_element_t<Matrix>;
200
201 if (col1 >= A.cols()) {
202 TH_MATH_ERROR("algebra::mat_swap_cols", col1, INVALID_ARGUMENT);
203 return mat_error(A);
204 }
205
206 if (col2 >= A.cols()) {
207 TH_MATH_ERROR("algebra::mat_swap_cols", col2, INVALID_ARGUMENT);
208 return mat_error(A);
209 }
210
211 if (col1 == col2)
212 return A;
213
214 for (unsigned int i = 0; i < A.cols(); ++i) {
215
216 const Type tmp = A(i, col1);
217
218 A(i, col1) = A(i, col2);
219 A(i, col2) = tmp;
220 }
221
222 return A;
223 }
224
225
232 template<typename Matrix, typename Type = matrix_element_t<Matrix>>
233 inline Matrix& mat_shift_diagonal(Matrix& A, const Type& sigma) {
234
235 const unsigned int count = min(A.rows(), A.cols());
236
237 for (unsigned int i = 0; i < count; ++i)
238 A(i, i) += sigma;
239
240 return A;
241 }
242
243
251 template<typename Type>
252 inline auto pair_inner_product(const Type& v_i, const Type& w_i) {
253 return v_i * w_i;
254 }
255
263 template<typename Type>
264 inline auto pair_inner_product(const complex<Type>& v_i, const complex<Type>& w_i) {
265 return v_i * conjugate(w_i);
266 }
267
268
274 template<typename Vector>
275 inline auto sqr_norm(const Vector& v) {
276
277 auto sum = pair_inner_product(v[0], v[0]);
278
279 // Use conjugation for complex numbers
280 for (unsigned int i = 1; i < v.size(); ++i)
281 sum += pair_inner_product(v[i], v[i]);
282
283 return sum;
284 }
285
286
291 template<typename Vector>
292 inline auto norm(const Vector& v) {
293
295 }
296
297
301 template<typename Vector>
302 inline Vector normalize(const Vector& v) {
303
304 Vector r;
305 r.resize(v.size());
306 vec_copy(r, v);
307
308 const auto m = norm(v);
309
310 if(abs(m) < MACH_EPSILON) {
311 TH_MATH_ERROR("algebra::normalize", m, DIV_BY_ZERO);
312 vec_error(r);
313 return r;
314 }
315
316 for (unsigned int i = 0; i < r.size(); ++i)
317 r[i] /= m;
318
319 return r;
320 }
321
322
326 template<typename Vector>
328
329 const auto m = norm(v);
330
331 if(abs(m) < MACH_EPSILON) {
332 TH_MATH_ERROR("algebra::make_normalized", m, DIV_BY_ZERO);
333 vec_error(v);
334 return v;
335 }
336
337 for (unsigned int i = 0; i < v.size(); ++i)
338 v[i] /= m;
339
340 return v;
341 }
342
343
350 template<typename Vector1, typename Vector2>
351 inline auto dot(const Vector1& v, const Vector2& w) {
352
353 if(v.size() != w.size()) {
354 TH_MATH_ERROR("algebra::dot", v.size(), INVALID_ARGUMENT);
356 }
357
358 auto sum = pair_inner_product(v[0], w[0]);
359
360 // Use conjugation for complex numbers
361 for (unsigned int i = 1; i < v.size(); ++i)
362 sum += pair_inner_product(v[i], w[i]);
363
364 return sum;
365 }
366
367
372 template<typename Vector1, typename Vector2>
373 inline Vector1 cross(const Vector1& v1, const Vector2& v2) {
374
375 Vector1 v3;
376 v3.resize(3);
377
378 if(v1.size() != 3) {
379 TH_MATH_ERROR("algebra::cross", v1.size(), INVALID_ARGUMENT);
380 vec_error(v3);
381 return v3;
382 }
383
384 if(v2.size() != 3) {
385 TH_MATH_ERROR("algebra::cross", v2.size(), INVALID_ARGUMENT);
386 vec_error(v3);
387 return v3;
388 }
389
390 v3[0] = v1[1] * v2[2] - v1[2] * v2[1];
391 v3[1] = v1[2] * v2[0] - v1[0] * v2[2];
392 v3[2] = v1[0] * v2[1] - v1[1] * v2[0];
393
394 return v3;
395 }
396
397
402 template<typename Matrix, typename MatrixT = Matrix>
403 inline MatrixT transpose(const Matrix& m) {
404
405 MatrixT res;
406 res.resize(m.cols(), m.rows());
407
408 for (unsigned int i = 0; i < m.rows(); ++i)
409 for (unsigned int j = 0; j < m.cols(); ++j)
410 res(j, i) = m(i, j);
411
412 return res;
413 }
414
415
419 template<typename Matrix>
421
422 if(m.rows() != m.cols()) {
423 TH_MATH_ERROR("algebra::make_transposed", m.rows(), INVALID_ARGUMENT);
424 mat_error(m);
425 return m;
426 }
427
428 for (unsigned int i = 0; i < m.rows(); ++i) {
429 for (unsigned int j = 0; j < i; ++j) {
430
431 const auto buff = m(i, j);
432 m(i, j) = m(j, i);
433 m(j, i) = buff;
434 }
435 }
436
437 return m;
438 }
439
440
447 template<typename Matrix1, typename Matrix2>
449
450 // Check that the two matrices have the correct
451 // number of rows and columns
452 if(src.rows() != dest.cols()) {
453 TH_MATH_ERROR("algebra::transpose", src.rows(), INVALID_ARGUMENT);
455 return dest;
456 }
457
458 if(src.cols() != dest.rows()) {
459 TH_MATH_ERROR("algebra::transpose", src.cols(), INVALID_ARGUMENT);
461 return dest;
462 }
463
464 // Overwrite dest with the transpose of src
465 for (unsigned int i = 0; i < src.rows(); ++i)
466 for (unsigned int j = 0; j < src.cols(); ++j)
467 dest(j, i) = src(i, j);
468
469 return dest;
470 }
471
472
477 template<typename Matrix, typename MatrixT = Matrix>
478 inline MatrixT hermitian(const Matrix& m) {
479
480 MatrixT res;
481 res.resize(m.cols(), m.rows());
482
483 for (unsigned int i = 0; i < m.rows(); ++i)
484 for (unsigned int j = 0; j < m.cols(); ++j)
485 res(j, i) = conjugate(m(i, j));
486
487 return res;
488 }
489
490
496 template<typename Matrix>
498
499 if(m.rows() != m.cols()) {
500 TH_MATH_ERROR("algebra::hermitian", m.rows(), INVALID_ARGUMENT);
501 mat_error(m);
502 return m;
503 }
504
505 for (unsigned int i = 0; i < m.rows(); ++i) {
506 for (unsigned int j = 0; j < i; ++j) {
507
508 const auto buff = m(i, j);
509 m(i, j) = conjugate(m(j, i));
510 m(j, i) = conjugate(buff);
511 }
512 }
513
514 return m;
515 }
516
517
525 template<typename Matrix1, typename Matrix2>
527
528 // Check that the two matrices have the correct
529 // number of rows and columns
530 if(src.rows() != dest.cols()) {
531 TH_MATH_ERROR("algebra::hermitian", src.rows(), INVALID_ARGUMENT);
533 return dest;
534 }
535
536 if(src.cols() != dest.rows()) {
537 TH_MATH_ERROR("algebra::hermitian", src.cols(), INVALID_ARGUMENT);
539 return dest;
540 }
541
542 // Overwrite dest with the transpose of src
543 for (unsigned int i = 0; i < src.rows(); ++i)
544 for (unsigned int j = 0; j < src.cols(); ++j)
545 dest(j, i) = conjugate(src(i, j));
546
547 return dest;
548 }
549
550
554 template<typename Matrix>
555 inline auto trace(const Matrix& m) {
556
557 auto sum = m(0, 0);
558 const size_t n = min(m.rows(), m.cols());
559
560 for (unsigned int i = 1; i < n; ++i)
561 sum += m(i, i);
562
563 return sum;
564 }
565
566
572 template<typename Matrix>
573 inline auto diagonal_product(const Matrix& m) {
574
575 auto mul = m(0, 0);
576 const size_t n = min(m.rows(), m.cols());
577
578 for (unsigned int i = 1; i < n; ++i)
579 mul *= m(i, i);
580
581 return mul;
582 }
583
584
589 template<typename Field, typename Matrix>
591
592 for (unsigned int i = 0; i < m.rows(); ++i)
593 for (unsigned int j = 0; j < m.cols(); ++j)
594 m(i, j) *= a;
595
596 return m;
597 }
598
599
606 template<typename Field, typename Matrix1, typename Matrix2>
608
609 if(src.rows() != dest.rows()) {
610 TH_MATH_ERROR("algebra::mat_scalmul", src.rows(), INVALID_ARGUMENT);
612 return dest;
613 }
614
615 if(src.cols() != dest.cols()) {
616 TH_MATH_ERROR("algebra::mat_scalmul", src.cols(), INVALID_ARGUMENT);
618 return dest;
619 }
620
621 for (unsigned int i = 0; i < src.rows(); ++i)
622 for (unsigned int j = 0; j < src.cols(); ++j)
623 dest(i, j) = a * src(i, j);
624
625 return dest;
626 }
627
628
633 template<typename Field, typename Vector>
635
636 for (unsigned int i = 0; i < v.size(); ++i)
637 v[i] *= a;
638
639 return v;
640 }
641
642
649 template<typename Field, typename Vector1, typename Vector2>
651
652 if(src.size() != dest.size()) {
653 TH_MATH_ERROR("algebra::vec_scalmul", src.size(), INVALID_ARGUMENT);
655 return dest;
656 }
657
658 for (unsigned int i = 0; i < src.size(); ++i)
659 dest[i] = a * src[i];
660
661 return dest;
662 }
663
664
665 // Operations involving a matrix and a vector
666
667
674 template<typename Matrix, typename Vector>
675 inline Vector& apply_transform(const Matrix& A, Vector& v) {
676
677 if(v.size() != A.cols()) {
678 TH_MATH_ERROR("algebra::apply_transform", v.size(), INVALID_ARGUMENT);
679 vec_error(v);
680 return v;
681 }
682
683 Vector res;
684 res.resize(v.size());
686
687 for (unsigned int i = 0; i < A.rows(); ++i)
688 for (unsigned int j = 0; j < A.cols(); ++j)
689 res[i] += A(i, j) * v[j];
690
691 vec_copy(res, v);
692 return v;
693 }
694
695
701 template<typename Matrix, typename Vector>
702 inline Vector transform(const Matrix& A, const Vector& v) {
703
704 Vector res;
705 res.resize(v.size());
706
707 if(v.size() != A.cols()) {
708 TH_MATH_ERROR("algebra::transform", v.size(), INVALID_ARGUMENT);
709 vec_error(res);
710 return res;
711 }
712
714
715 for (unsigned int i = 0; i < A.rows(); ++i)
716 for (unsigned int j = 0; j < A.cols(); ++j)
717 res[i] += A(i, j) * v[j];
718
719 return res;
720 }
721
722
730 template<typename Matrix, typename Vector1, typename Vector2>
731 inline Vector1& transform(Vector1& res, const Matrix& A, const Vector2& v) {
732
733 if(v.size() != A.cols()) {
734 TH_MATH_ERROR("algebra::transform", v.size(), INVALID_ARGUMENT);
735 vec_error(res);
736 return res;
737 }
738
739 if(res.size() != v.size()) {
740 TH_MATH_ERROR("algebra::transform", res.size(), INVALID_ARGUMENT);
741 vec_error(res);
742 return res;
743 }
744
746
747 for (unsigned int i = 0; i < A.rows(); ++i)
748 for (unsigned int j = 0; j < A.cols(); ++j)
749 res[i] += A(i, j) * v[j];
750
751 return res;
752 }
753
754
755 // Operations involving multiple matrices or vectors
756
757
763 template<typename Matrix1, typename Matrix2>
764 inline Matrix2& mat_sum(Matrix1& A, const Matrix2& B) {
765
766 if(A.rows() != B.rows()) {
767 TH_MATH_ERROR("algebra::mat_sum", A.rows(), INVALID_ARGUMENT);
768 mat_error(A);
769 return A;
770 }
771
772 if(A.cols() != B.cols()) {
773 TH_MATH_ERROR("algebra::mat_sum", A.cols(), INVALID_ARGUMENT);
774 mat_error(A);
775 return A;
776 }
777
778 for (unsigned int i = 0; i < A.rows(); ++i)
779 for (unsigned int j = 0; j < A.cols(); ++j)
780 A(i, j) = A(i, j) + B(i, j);
781
782 return A;
783 }
784
785
791 template<typename Matrix1, typename Matrix2, typename Matrix3>
792 inline Matrix1& mat_sum(Matrix1& res, const Matrix2& A, const Matrix3& B) {
793
794 if(A.rows() != B.rows()) {
795 TH_MATH_ERROR("algebra::mat_sum", A.rows(), INVALID_ARGUMENT);
796 mat_error(res);
797 return res;
798 }
799
800 if(A.cols() != B.cols()) {
801 TH_MATH_ERROR("algebra::mat_sum", A.cols(), INVALID_ARGUMENT);
802 mat_error(res);
803 return res;
804 }
805
806 if(res.rows() != A.rows()) {
807 TH_MATH_ERROR("algebra::mat_sum", res.rows(), INVALID_ARGUMENT);
808 mat_error(res);
809 return res;
810 }
811
812 if(res.cols() != A.cols()) {
813 TH_MATH_ERROR("algebra::mat_sum", res.cols(), INVALID_ARGUMENT);
814 mat_error(res);
815 return res;
816 }
817
818 for (unsigned int i = 0; i < A.rows(); ++i)
819 for (unsigned int j = 0; j < A.cols(); ++j)
820 res(i, j) = A(i, j) + B(i, j);
821
822 return res;
823 }
824
825
831 template<typename Matrix1, typename Matrix2>
832 inline Matrix2& mat_diff(Matrix1& A, const Matrix2& B) {
833
834 if(A.rows() != B.rows()) {
835 TH_MATH_ERROR("algebra::mat_diff", A.rows(), INVALID_ARGUMENT);
836 mat_error(A);
837 return A;
838 }
839
840 if(A.cols() != B.cols()) {
841 TH_MATH_ERROR("algebra::mat_diff", A.cols(), INVALID_ARGUMENT);
842 mat_error(A);
843 return A;
844 }
845
846 for (unsigned int i = 0; i < A.rows(); ++i)
847 for (unsigned int j = 0; j < A.cols(); ++j)
848 A(i, j) = A(i, j) - B(i, j);
849
850 return A;
851 }
852
853
859 template<typename Matrix1, typename Matrix2, typename Matrix3>
860 inline Matrix1& mat_diff(Matrix1& res, const Matrix2& A, const Matrix3& B) {
861
862 if(A.rows() != B.rows()) {
863 TH_MATH_ERROR("algebra::mat_diff", A.rows(), INVALID_ARGUMENT);
864 mat_error(res);
865 return res;
866 }
867
868 if(A.cols() != B.cols()) {
869 TH_MATH_ERROR("algebra::mat_diff", A.cols(), INVALID_ARGUMENT);
870 mat_error(res);
871 return res;
872 }
873
874 if(res.rows() != A.rows()) {
875 TH_MATH_ERROR("algebra::mat_diff", res.rows(), INVALID_ARGUMENT);
876 mat_error(res);
877 return res;
878 }
879
880 if(res.cols() != A.cols()) {
881 TH_MATH_ERROR("algebra::mat_diff", res.cols(), INVALID_ARGUMENT);
882 mat_error(res);
883 return res;
884 }
885
886 for (unsigned int i = 0; i < A.rows(); ++i)
887 for (unsigned int j = 0; j < A.cols(); ++j)
888 res(i, j) = A(i, j) - B(i, j);
889
890 return res;
891 }
892
893
900 template<typename Field1, typename Matrix1, typename Field2, typename Matrix2>
902 Field1 alpha, Matrix1& A, Field2 beta, const Matrix2& B) {
903
904 if(A.rows() != B.rows()) {
905 TH_MATH_ERROR("algebra::mat_lincomb", A.rows(), INVALID_ARGUMENT);
906 mat_error(A);
907 return A;
908 }
909
910 if(A.cols() != B.cols()) {
911 TH_MATH_ERROR("algebra::mat_lincomb", A.cols(), INVALID_ARGUMENT);
912 mat_error(A);
913 return A;
914 }
915
916 for (unsigned int i = 0; i < A.rows(); ++i)
917 for (unsigned int j = 0; j < A.cols(); ++j)
918 A(i, j) = A(i, j) * alpha + B(i, j) * beta;
919
920 return A;
921 }
922
923
933 template<typename Matrix1, typename Field1, typename Matrix2,
934 typename Field2, typename Matrix3>
936 Matrix1& res, Field1 alpha, const Matrix2& A, Field2 beta, const Matrix3& B) {
937
938 if(A.rows() != B.rows()) {
939 TH_MATH_ERROR("algebra::mat_lincomb", A.rows(), INVALID_ARGUMENT);
940 mat_error(A);
941 return A;
942 }
943
944 if(A.cols() != B.cols()) {
945 TH_MATH_ERROR("algebra::mat_lincomb", A.cols(), INVALID_ARGUMENT);
946 mat_error(A);
947 return A;
948 }
949
950 if(res.rows() != A.rows()) {
951 TH_MATH_ERROR("algebra::mat_lincomb", res.rows(), INVALID_ARGUMENT);
952 mat_error(res);
953 return res;
954 }
955
956 if(res.cols() != A.cols()) {
957 TH_MATH_ERROR("algebra::mat_lincomb", res.cols(), INVALID_ARGUMENT);
958 mat_error(res);
959 return res;
960 }
961
962 for (unsigned int i = 0; i < A.rows(); ++i)
963 for (unsigned int j = 0; j < A.cols(); ++j)
964 res(i, j) = A(i, j) * alpha + B(i, j) * beta;
965
966 return A;
967 }
968
969
976 template<typename Matrix1, typename Matrix2, typename Matrix3 = Matrix1>
977 inline Matrix3 mat_mul(const Matrix1& A, const Matrix2& B) {
978
979 Matrix3 R;
980 R.resize(A.rows(), B.cols());
981
982 if(A.cols() != B.rows()) {
983 TH_MATH_ERROR("algebra::mat_mul", A.cols(), INVALID_ARGUMENT);
984 mat_error(R);
985 return R;
986 }
987
988 mat_zeroes(R);
989
990 for (unsigned int i = 0; i < A.rows(); ++i)
991 for (unsigned int j = 0; j < B.cols(); ++j)
992 for (unsigned int k = 0; k < A.cols(); ++k)
993 R(i, j) += A(i, k) * B(k, j);
994
995 return R;
996 }
997
998
1006 template<typename Matrix1, typename Matrix2, typename Matrix3>
1007 inline Matrix1& mat_mul(Matrix1& R, const Matrix2& A, const Matrix3& B) {
1008
1009 if(R.rows() != A.rows()) {
1010 TH_MATH_ERROR("algebra::mat_mul", R.rows(), INVALID_ARGUMENT);
1011 mat_error(R);
1012 return R;
1013 }
1014
1015 if(R.cols() != B.cols()) {
1016 TH_MATH_ERROR("algebra::mat_mul", R.cols(), INVALID_ARGUMENT);
1017 mat_error(R);
1018 return R;
1019 }
1020
1021 if(A.cols() != B.rows()) {
1022 TH_MATH_ERROR("algebra::mat_mul", A.cols(), INVALID_ARGUMENT);
1023 mat_error(R);
1024 return R;
1025 }
1026
1027 mat_zeroes(R);
1028
1029 for (unsigned int i = 0; i < A.rows(); ++i)
1030 for (unsigned int j = 0; j < B.cols(); ++j)
1031 for (unsigned int k = 0; k < A.cols(); ++k)
1032 R(i, j) += A(i, k) * B(k, j);
1033
1034 return R;
1035 }
1036
1037
1048 template<typename Matrix1, typename Matrix2, typename Matrix3 = Matrix1>
1049 inline Matrix3 mat_transpose_mul(const Matrix1& A, const Matrix2& B) {
1050
1051 Matrix3 R;
1052 R.resize(A.cols(), B.cols());
1053
1054 if(A.rows() != B.rows()) {
1055 TH_MATH_ERROR("algebra::mat_transpose_mul", A.rows(), INVALID_ARGUMENT);
1056 return mat_error(R);
1057 }
1058
1059 mat_zeroes(R);
1060
1061 for (unsigned int i = 0; i < R.rows(); ++i)
1062 for (unsigned int j = 0; j < R.cols(); ++j)
1063 for (unsigned int k = 0; k < A.rows(); ++k)
1064 R(i, j) += A(k, i) * B(k, j);
1065
1066 return R;
1067 }
1068
1069
1080 template<typename Matrix1, typename Matrix2, typename Matrix3 = Matrix1>
1081 inline Matrix3 mat_mul_transpose(const Matrix1& A, const Matrix2& B) {
1082
1083 Matrix3 R;
1084 R.resize(A.rows(), B.rows());
1085
1086 if(A.cols() != B.cols()) {
1087 TH_MATH_ERROR("algebra::mat_mul_transpose", A.rows(), INVALID_ARGUMENT);
1088 return mat_error(R);
1089 }
1090
1091 mat_zeroes(R);
1092
1093 for (unsigned int i = 0; i < R.rows(); ++i)
1094 for (unsigned int j = 0; j < R.cols(); ++j)
1095 for (unsigned int k = 0; k < A.cols(); ++k)
1096 R(i, j) += A(i, k) * B(j, k);
1097
1098 return R;
1099 }
1100
1101
1108 template<typename Matrix1, typename Matrix2>
1109 inline bool mat_equals(
1110 const Matrix1& A, const Matrix2& B, real tolerance = 10 * MACH_EPSILON) {
1111
1112 if(A.rows() != B.rows() || A.cols() != B.cols())
1113 return false;
1114
1115 for (unsigned int i = 0; i < A.rows(); ++i) {
1116 for (unsigned int j = 0; j < A.cols(); ++j) {
1117
1118 const auto diff = abs(A(i, j) - B(i, j));
1119
1120 if(diff > tolerance || is_nan(diff))
1121 return false;
1122 }
1123 }
1124
1125 return true;
1126 }
1127
1128
1134 template<typename Vector1, typename Vector2>
1135 inline Vector2& vec_sum(Vector1& v1, const Vector2& v2) {
1136
1137 if(v1.size() != v2.size()) {
1138 TH_MATH_ERROR("algebra::vec_sum", v1.size(), INVALID_ARGUMENT);
1139 vec_error(v1);
1140 return v1;
1141 }
1142
1143 for (unsigned int i = 0; i < v1.size(); ++i)
1144 v1[i] = v1[i] + v2[i];
1145
1146 return v1;
1147 }
1148
1149
1155 template<typename Vector1, typename Vector2, typename Vector3>
1156 inline Vector1& vec_sum(Vector1& res, const Vector2& v1, const Vector3& v2) {
1157
1158 if(v1.size() != v2.size()) {
1159 TH_MATH_ERROR("algebra::vec_sum", v1.size(), INVALID_ARGUMENT);
1160 vec_error(res);
1161 return res;
1162 }
1163
1164 if(res.size() != v1.size()) {
1165 TH_MATH_ERROR("algebra::vec_sum", res.size(), INVALID_ARGUMENT);
1166 vec_error(res);
1167 return res;
1168 }
1169
1170 for (unsigned int i = 0; i < v1.size(); ++i)
1171 res[i] = v1[i] + v2[i];
1172
1173 return res;
1174 }
1175
1176
1182 template<typename Vector1, typename Vector2>
1183 inline Vector2& vec_diff(Vector1& v1, const Vector2& v2) {
1184
1185 if(v1.size() != v2.size()) {
1186 TH_MATH_ERROR("algebra::vec_diff", v1.size(), INVALID_ARGUMENT);
1187 vec_error(v1);
1188 return v1;
1189 }
1190
1191 for (unsigned int i = 0; i < v1.size(); ++i)
1192 v1[i] = v1[i] - v2[i];
1193
1194 return v1;
1195 }
1196
1197
1203 template<typename Vector1, typename Vector2, typename Vector3>
1204 inline Vector1& vec_diff(Vector1& res, const Vector2& v1, const Vector3& v2) {
1205
1206 if(v1.size() != v2.size()) {
1207 TH_MATH_ERROR("algebra::vec_diff", v1.size(), INVALID_ARGUMENT);
1208 vec_error(res);
1209 return res;
1210 }
1211
1212 if(res.size() != v1.size()) {
1213 TH_MATH_ERROR("algebra::vec_diff", res.size(), INVALID_ARGUMENT);
1214 vec_error(res);
1215 return res;
1216 }
1217
1218 for (unsigned int i = 0; i < v1.size(); ++i)
1219 res[i] = v1[i] - v2[i];
1220
1221 return res;
1222 }
1223
1224
1225 // Matrix properties
1226
1227
1231 template<typename Matrix>
1232 inline bool is_square(const Matrix& m) {
1233 return (m.rows() == m.cols());
1234 }
1235
1236
1242 template<typename Matrix>
1243 inline bool is_diagonal(const Matrix& m, real tolerance = 10 * MACH_EPSILON) {
1244
1245 for (unsigned int i = 0; i < m.rows(); ++i)
1246 for (unsigned int j = 0; j < m.cols(); ++j)
1247 if(i != j && abs(m(i, j)) > tolerance)
1248 return false;
1249
1250 return true;
1251 }
1252
1258 template<typename Matrix>
1260
1261 if(!is_square(m))
1262 return false;
1263
1264 for (unsigned int i = 0; i < m.rows(); ++i)
1265 for (unsigned int j = 0; j < m.cols(); ++j)
1266 if(abs(m(i, j) - m(j, i)) > tolerance)
1267 return false;
1268
1269 return true;
1270 }
1271
1272
1278 template<typename Matrix>
1280
1281 if(!is_square(m))
1282 return false;
1283
1284 for (unsigned int i = 0; i < m.rows(); ++i)
1285 for (unsigned int j = i + 1; j < m.cols(); ++j)
1286 if (abs(m(i, j)) > tolerance)
1287 return false;
1288
1289 return true;
1290 }
1291
1292
1298 template<typename Matrix>
1300
1301 if(!is_square(m))
1302 return false;
1303
1304 for (unsigned int i = 0; i < m.rows(); ++i)
1305 for (unsigned int j = 0; j < i; ++j)
1306 if (abs(m(i, j)) > tolerance)
1307 return false;
1308
1309 return true;
1310 }
1311
1312
1322 template<typename Matrix, enable_matrix<Matrix> = true>
1323 inline real density(const Matrix& A, real tolerance = 1E-12) {
1324
1325 long unsigned int n = 0;
1326
1327 for (unsigned int i = 0; i < A.rows(); ++i)
1328 for (unsigned int j = 0; j < A.cols(); ++j)
1329 if (abs(A(i, j)) > tolerance)
1330 n++;
1331
1332 return (real(n) / A.rows()) / A.cols();
1333 }
1334
1335
1345 template<typename Matrix, enable_matrix<Matrix> = true>
1346 inline real sparsity(const Matrix& A, real tolerance = 1E-12) {
1347
1348 return 1.0 - density(A, tolerance);
1349 }
1350
1351
1352 // Matrix decompositions
1353
1354
1362 template<typename Matrix1, typename Matrix2, typename Matrix3>
1363 inline void decompose_lu(const Matrix1& A, Matrix2& L, Matrix3& U) {
1364
1365 // Check the shapes of A, L and U
1366 unsigned int err = 0;
1367
1368 // Make sure to allocate L and U if they are empty
1369 L.resize(A.rows(), A.cols());
1370 U.resize(A.rows(), A.cols());
1371
1372 if (!is_square(A))
1373 err = A.rows();
1374
1375 if (A.rows() != L.rows())
1376 err = L.rows();
1377
1378 if (A.cols() != L.cols())
1379 err = L.cols();
1380
1381 if (A.rows() != U.rows())
1382 err = U.rows();
1383
1384 if (A.cols() != U.cols())
1385 err = U.cols();
1386
1387 if (err) {
1388 TH_MATH_ERROR("algebra::decompose_lu", err, INVALID_ARGUMENT);
1390 return;
1391 }
1392
1394
1395 // Set the diagonal of L to 1.0
1396 for (unsigned int i = 0; i < A.rows(); ++i)
1397 L(i, i) = (Type) 1.0;
1398
1399 // Compute L and U
1400 for(unsigned int i = 0; i < A.rows(); ++i) {
1401
1402 for(unsigned int j = 0; j < A.rows(); ++j) {
1403
1404 U(i, j) = A(i, j);
1405
1406 for(unsigned int k = 0; k < i; ++k)
1407 U(i, j) -= pair_inner_product(L(i, k), U(k, j));
1408 }
1409
1410 for(unsigned int j = i + 1; j < A.rows(); ++j) {
1411
1412 L(j, i) = A(j, i);
1413
1414 for(unsigned int k = 0; k < i; ++k)
1415 L(j, i) -= pair_inner_product(L(j, k), U(k, i));
1416
1417 L(j, i) /= U(i, i);
1418 }
1419 }
1420 }
1421
1422
1431 template<typename Matrix>
1433
1434 if (!is_square(A)) {
1435 TH_MATH_ERROR("algebra::decompose_lu_inplace", A.rows(), INVALID_ARGUMENT);
1436 mat_error(A);
1437 return A;
1438 }
1439
1440 for(unsigned int j = 0; j < A.cols(); ++j) {
1441
1442 for(unsigned int i = j + 1; i < A.rows(); ++i) {
1443
1444 A(i, j) /= A(j, j);
1445
1446 for(unsigned int k = j + 1; k < A.rows(); ++k)
1447 A(i, k) -= pair_inner_product(A(i, j), A(j, k));
1448 }
1449 }
1450
1451 return A;
1452 }
1453
1454
1461 template<typename Matrix>
1463
1464 Matrix L;
1465 L.resize(A.rows(), A.cols());
1467
1468 if (!is_square(A)) {
1469 TH_MATH_ERROR("algebra::decompose_cholesky", A.rows(), INVALID_ARGUMENT);
1470 return mat_error(L);
1471 }
1472
1473 if (!is_symmetric(A)) {
1474 TH_MATH_ERROR("algebra::decompose_cholesky", false, INVALID_ARGUMENT);
1475 return mat_error(L);
1476 }
1477
1478 mat_zeroes(L);
1479
1480 for (unsigned int i = 0; i < A.cols(); ++i) {
1481
1482 for (unsigned int j = 0; j <= i; ++j) {
1483
1484 Type sum = pair_inner_product(L(i, 0), L(j, 0));
1485
1486 for (unsigned int k = 1; k < j; ++k)
1487 sum += pair_inner_product(L(i, k), L(j, k));
1488
1489 if (i == j) {
1490
1491 const Type sqr_diag = A(j, j) - sum;
1492
1493 // Additional check to ensure that the matrix is positive definite
1494 if (sqr_diag < MACH_EPSILON) {
1495 TH_MATH_ERROR("algebra::decompose_cholesky", sqr_diag, INVALID_ARGUMENT);
1496 return mat_error(L);
1497 }
1498
1499 L(i, j) = sqrt(sqr_diag);
1500
1501 } else {
1502
1503 L(i, j) = (A(i, j) - sum) / L(j, j);
1504 }
1505 }
1506 }
1507
1508 return L;
1509 }
1510
1511
1517 template<typename Matrix>
1519
1521
1522 if (!is_square(A)) {
1523 TH_MATH_ERROR("algebra::decompose_cholesky_inplace", A.rows(), INVALID_ARGUMENT);
1524 return mat_error(A);
1525 }
1526
1527 if (!is_symmetric(A)) {
1528 TH_MATH_ERROR("algebra::decompose_cholesky_inplace", false, INVALID_ARGUMENT);
1529 return mat_error(A);
1530 }
1531
1532 // Compute the Cholesky decomposition in-place
1533 for (unsigned int k = 0; k < A.rows(); ++k) {
1534
1535 A(k, k) = sqrt(A(k, k));
1536
1537 for (unsigned int i = k + 1; i < A.rows(); ++i)
1538 A(i, k) = A(i, k) / A(k, k);
1539
1540 for (unsigned int j = k + 1; j < A.rows(); ++j)
1541 for (unsigned int i = j; i < A.cols(); ++i)
1542 A(i, j) -= pair_inner_product(A(i, k), A(j, k));
1543 }
1544
1545 // Zero out elements over the diagonal
1546 for (unsigned int i = 0; i < A.rows(); ++i)
1547 for (unsigned int j = i + 1; j < A.rows(); ++j)
1548 A(i, j) = Type(0.0);
1549
1550
1551 return A;
1552 }
1553
1554
1555 // Linear system solvers
1556
1557
1563 template<typename Matrix, typename Vector>
1564 inline Vector solve_triangular_lower(const Matrix& L, const Vector& b) {
1565
1566 Vector x;
1567 x.resize(L.cols());
1569
1570 if (!is_square(L)) {
1571 TH_MATH_ERROR("algebra::solve_triangular_lower", false, INVALID_ARGUMENT);
1572 return vec_error(x);
1573 }
1574
1575 if (b.size() != L.rows()) {
1576 TH_MATH_ERROR("algebra::solve_triangular_lower", b.size(), INVALID_ARGUMENT);
1577 return vec_error(x);
1578 }
1579
1580 // Solve using forward substitution
1581 for (unsigned int i = 0; i < L.cols(); ++i) {
1582
1583 Type sum = L(i, 0) * x[0];
1584
1585 for (unsigned int j = 1; j < i; ++j)
1586 sum += L(i, j) * x[j];
1587
1588 if (abs(L(i, i)) < MACH_EPSILON) {
1589 TH_MATH_ERROR("algebra::solve_triangular_lower", L(i, i), DIV_BY_ZERO);
1590 return vec_error(x);
1591 }
1592
1593 x[i] = (b[i] - sum) / L(i, i);
1594 }
1595
1596 return x;
1597 }
1598
1599
1605 template<typename Matrix, typename Vector>
1606 inline Vector solve_triangular_upper(const Matrix& U, const Vector& b) {
1607
1608 Vector x;
1609 x.resize(U.cols());
1611
1612 if (!is_square(U)) {
1613 TH_MATH_ERROR("solve_triangular_upper", false, INVALID_ARGUMENT);
1614 return vec_error(x);
1615 }
1616
1617 if (b.size() != U.rows()) {
1618 TH_MATH_ERROR("solve_triangular_upper", b.size(), INVALID_ARGUMENT);
1619 return vec_error(x);
1620 }
1621
1622 if (U.rows() == 1) {
1623 x[0] = b[0] / U(0, 0);
1624 return x;
1625 }
1626
1627 // Solve using backward substitution
1628 for (int i = U.rows() - 1; i >= 0; --i) {
1629
1630 Type sum = U(i, i + 1) * x[i + 1];
1631
1632 for (unsigned int j = i + 2; j < U.cols(); ++j)
1633 sum += U(i, j) * x[j];
1634
1635 if (abs(U(i, i)) < MACH_EPSILON) {
1636 TH_MATH_ERROR("solve_triangular_upper", U(i, i), DIV_BY_ZERO);
1637 return vec_error(x);
1638 }
1639
1640 x[i] = (b[i] - sum) / U(i, i);
1641 }
1642
1643 return x;
1644 }
1645
1646
1654 template<typename Matrix, typename Vector>
1655 inline Vector solve_triangular(const Matrix& T, const Vector& b) {
1656
1657 // Pick the correct solver
1659 return solve_triangular_lower(T, b);
1660 else if(is_upper_triangular(T))
1661 return solve_triangular_upper(T, b);
1662 else {
1663 Vector err;
1664 err.resize(b.size());
1665 return vec_error(err);
1666 }
1667 }
1668
1669
1680 template<typename Matrix, typename Vector>
1681 inline Vector& solve_lu_inplace(const Matrix& A, Vector& b) {
1682
1683 if (!is_square(A)) {
1684 TH_MATH_ERROR("algebra::solve_lu_inplace", A.rows(), INVALID_ARGUMENT);
1685 return vec_error(b);
1686 }
1687
1688 if (A.rows() != b.size()) {
1689 TH_MATH_ERROR("algebra::solve_lu_inplace", A.rows(), INVALID_ARGUMENT);
1690 return vec_error(b);
1691 }
1692
1694
1695 // Forward elimination: solves the lower system (L matrix).
1696 for (unsigned int i = 1; i < A.rows(); ++i) {
1697
1698 Type sum = Type(0.0);
1699
1700 for (unsigned int j = 0; j < i; ++j)
1701 sum += A(i, j) * b[j];
1702
1703 b[i] -= sum;
1704 }
1705
1706 // Backward elimination: solves the upper system (U matrix).
1707 for (int i = A.rows() - 1; i >= 0; --i) {
1708
1709 Type sum = Type(0.0);
1710
1711 for (unsigned int j = i + 1; j < A.rows(); ++j)
1712 sum += A(i, j) * b[j];
1713
1714 b[i] = (b[i] - sum) / A(i, i);
1715 }
1716
1717 return b;
1718 }
1719
1720
1728 template<typename Matrix, typename Vector>
1730
1731 // Apply in-place LU decomposition
1733
1734 // Apply forward and backward substitution
1735 return solve_lu_inplace(A, b);
1736 }
1737
1738
1749 template<typename Matrix1, typename Matrix2, typename Vector>
1750 inline Vector solve_lu(const Matrix1& L, const Matrix2& U, const Vector& b) {
1751
1752 Vector x = b;
1753
1754 if (!is_square(L)) {
1755 TH_MATH_ERROR("algebra::solve_lu", L.rows(), INVALID_ARGUMENT);
1756 return vec_error(x);
1757 }
1758
1759 if (!is_square(U)) {
1760 TH_MATH_ERROR("algebra::solve_lu", U.rows(), INVALID_ARGUMENT);
1761 return vec_error(x);
1762 }
1763
1764 if (L.rows() != U.rows()) {
1765 TH_MATH_ERROR("algebra::solve_lu", U.rows(), INVALID_ARGUMENT);
1766 return vec_error(x);
1767 }
1768
1769 if (b.size() != L.rows()) {
1770 TH_MATH_ERROR("algebra::solve_lu", b.size(), INVALID_ARGUMENT);
1771 return vec_error(x);
1772 }
1773
1775
1776 // Forward elimination for L
1777 for (unsigned int i = 1; i < L.rows(); ++i) {
1778
1779 Type1 sum = Type1(0.0);
1780
1781 for (unsigned int j = 0; j < i; ++j)
1782 sum += L(i, j) * x[j];
1783
1784 x[i] -= sum;
1785 }
1786
1788
1789 // Backward elimination for U
1790 for (int i = U.rows() - 1; i >= 0; --i) {
1791
1792 Type2 sum = Type2(0.0);
1793
1794 for (unsigned int j = i + 1; j < U.rows(); ++j)
1795 sum += U(i, j) * x[j];
1796
1797 if (abs(U(i, i)) < MACH_EPSILON) {
1798 TH_MATH_ERROR("algebra::solve_lu", U(i, i), DIV_BY_ZERO);
1799 return vec_error(x);
1800 }
1801
1802 x[i] = (x[i] - sum) / U(i, i);
1803 }
1804
1805 return x;
1806 }
1807
1808
1817 template<typename Matrix, typename Vector>
1818 inline Vector solve_cholesky(const Matrix& L, const Vector& b) {
1819
1820 Vector x = b;
1821
1822 if (!is_square(L)) {
1823 TH_MATH_ERROR("algebra::solve_cholesky", L.rows(), INVALID_ARGUMENT);
1824 return vec_error(x);
1825 }
1826
1827 if (L.rows() != b.size()) {
1828 TH_MATH_ERROR("algebra::solve_cholesky", b.size(), INVALID_ARGUMENT);
1829 return vec_error(x);
1830 }
1831
1833
1834 // Forward elimination for L
1835 for (unsigned int i = 0; i < L.rows(); ++i) {
1836
1837 Type sum = Type(0.0);
1838
1839 for (unsigned int j = 0; j < i; ++j)
1840 sum += L(i, j) * x[j];
1841
1842 x[i] = (x[i] - sum) / L(i, i);
1843 }
1844
1845 // Backward elimination for L transpose
1846 for (int i = L.rows() - 1; i >= 0; --i) {
1847
1848 Type sum = Type(0.0);
1849
1850 for (unsigned int j = i + 1; j < L.rows(); ++j)
1851 sum += L(j, i) * x[j];
1852
1853 x[i] = (x[i] - sum) / L(i, i);
1854 }
1855
1856 return x;
1857 }
1858
1859
1866 template<typename Matrix, typename Vector>
1867 inline Vector solve(const Matrix& A, const Vector& b) {
1868
1869 // Use LU decomposition to solve the linear system
1870 return solve_lu(A, b);
1871 }
1872
1873
1874 // Other composite operations
1875
1876
1882 template<typename Matrix1, typename Matrix2>
1884
1885 if(src.rows() != src.cols()) {
1886 TH_MATH_ERROR("algebra::inverse", src.rows(), INVALID_ARGUMENT);
1887 mat_error(dest);
1888 return dest;
1889 }
1890
1891 if(dest.rows() != src.rows()) {
1892 TH_MATH_ERROR("algebra::inverse", dest.rows(), INVALID_ARGUMENT);
1893 mat_error(dest);
1894 return dest;
1895 }
1896
1897 if(dest.cols() != src.cols()) {
1898 TH_MATH_ERROR("algebra::inverse", dest.cols(), INVALID_ARGUMENT);
1899 mat_error(dest);
1900 return dest;
1901 }
1902
1904
1905 // Prepare extended matrix (A|B)
1906 Matrix1 A;
1907 A.resize(src.rows(), src.cols());
1908 dest.resize(src.rows(), src.cols());
1909 mat_copy(A, src);
1911
1912 // Iterate on all columns
1913 for (unsigned int i = 0; i < src.rows(); ++i) {
1914
1915 // Make sure the element on the diagonal
1916 // is non-zero by adding the first non-zero row
1917 if(A(i, i) == (Type) 0) {
1918
1919 bool flag = false;
1920
1921 // Iterate on all rows
1922 for (unsigned int j = i + 1; j < src.rows(); ++j) {
1923
1924 // Add the j-th row to the i-th row
1925 // if Aji is non-zero
1926 if(A(j, i) != (Type) 0) {
1927
1928 for (unsigned int k = 0; k < src.rows(); ++k) {
1929 A(i, k) += A(j, k);
1930 dest(i, k) += dest(j, k);
1931 }
1932
1933 flag = true;
1934 break;
1935 }
1936 }
1937
1938 if(!flag) {
1939 TH_MATH_ERROR("algebra::inverse", flag, IMPOSSIBLE_OPERATION);
1940 mat_error(dest);
1941 return dest;
1942 }
1943 }
1944
1945 auto inv_pivot = ((Type) 1.0) / A(i, i);
1946
1947 // Use the current row to make all other
1948 // elements of the column equal to zero
1949 for (unsigned int j = 0; j < src.rows(); ++j) {
1950
1951 // Skip the current row
1952 if(j == i)
1953 continue;
1954
1955 // Multiplication coefficient for
1956 // the elision of Ajk
1957 const auto coeff = A(j, i) * inv_pivot;
1958
1959 for (unsigned int k = 0; k < src.rows(); ++k) {
1960 A(j, k) -= coeff * A(i, k);
1961 dest(j, k) -= coeff * dest(i, k);
1962 }
1963 }
1964
1965 // Divide the current row by the pivot
1966 for (unsigned int j = 0; j < src.rows(); ++j) {
1967 A(i, j) *= inv_pivot;
1968 dest(i, j) *= inv_pivot;
1969 }
1970
1971 }
1972
1973 return dest;
1974 }
1975
1976
1982 template<typename Matrix, typename MatrixInv = Matrix>
1983 inline MatrixInv inverse(const Matrix& m) {
1984 MatrixInv res;
1985 res.resize(m.rows(), m.cols());
1986 inverse(res, m);
1987 return res;
1988 }
1989
1990
1996 template<typename Matrix>
1997 inline Matrix& invert(Matrix& m) {
1998
1999 if(m.rows() != m.cols()) {
2000 TH_MATH_ERROR("algebra::invert", m.rows(), INVALID_ARGUMENT);
2001 mat_error(m);
2002 return m;
2003 }
2004
2005 // Prepare extended matrix (A|B)
2006 Matrix tmp;
2007 tmp.resize(m.rows(), m.cols());
2008 inverse(tmp, m);
2009
2010 // Modify the matrix only when the inversion
2011 // has succeeded
2012 mat_copy(m, tmp);
2013 return m;
2014 }
2015
2016
2023 template<typename Matrix>
2024 inline auto det(const Matrix& A) {
2025
2026 Matrix LU = A;
2028
2029 // The determinant of a triangular matrix
2030 // is the product of the elements on its diagonal
2031 return diagonal_product(LU);
2032 }
2033
2034
2035 // Eigensolvers
2036
2037
2045 template<typename Matrix, typename Vector>
2046 inline auto rayleigh_quotient(const Matrix& A, const Vector& x) {
2047
2048 const auto p = dot(x, x);
2049
2050 // Check for division by zero
2051 if (abs(p) < MACH_EPSILON) {
2052 TH_MATH_ERROR("rayleigh_quotient", abs(p), DIV_BY_ZERO);
2053 return vector_element_t<Vector>(nan());
2054 }
2055
2056 return dot(x, transform(A, x)) / p;
2057 }
2058
2059
2071 template<typename Matrix, typename Vector>
2072 inline auto eigenvalue_power(
2073 const Matrix& A, const Vector& x,
2075
2077
2078 if (!is_square(A)) {
2079 TH_MATH_ERROR("eigenvalue_power", is_square(A), INVALID_ARGUMENT);
2080 return Type(nan());
2081 }
2082
2083 if (x.size() != A.rows()) {
2084 TH_MATH_ERROR("eigenvalue_power", x.size(), INVALID_ARGUMENT);
2085 return Type(nan());
2086 }
2087
2088 // Apply a first iteration to initialize
2089 // the current and previous vectors
2090 Vector x_prev = x;
2092
2093 // Iteration counter
2094 unsigned int i;
2095
2096 // Iteratively apply the matrix to the vector
2097 for (i = 1; i <= max_iter; ++i) {
2098
2099 x_prev = x_curr;
2101
2102 // Stop the algorithm when |x_k+1 +- x_k| is
2103 // less then the tolerance in module
2104 if (norm(x_curr - x_prev) <= tolerance ||
2106 break;
2107 }
2108
2109 // The algorithm did not converge
2110 if (i > max_iter) {
2111 TH_MATH_ERROR("eigenvalue_power", i, NO_ALGO_CONVERGENCE);
2112 return Type(nan());
2113 }
2114
2115 return dot(x_curr, A * x_curr);
2116 }
2117
2118
2132 template<typename Matrix, typename Vector1, typename Vector2 = Vector1>
2133 inline auto eigenpair_power(
2134 const Matrix& A, const Vector1& x, Vector2& v,
2136
2138
2139 if (!is_square(A)) {
2140 TH_MATH_ERROR("eigenpair_power", is_square(A), INVALID_ARGUMENT);
2141 return Type(nan());
2142 }
2143
2144 if (x.size() != A.rows()) {
2145 TH_MATH_ERROR("eigenpair_power", x.size(), INVALID_ARGUMENT);
2146 return Type(nan());
2147 }
2148
2149 if (v.size() != x.size()) {
2150 TH_MATH_ERROR("eigenpair_power", v.size(), INVALID_ARGUMENT);
2151 return Type(nan());
2152 }
2153
2154 // Apply a first iteration to initialize
2155 // the current and previous vectors
2156 Vector1 x_prev = x;
2158
2159 // Iteration counter
2160 unsigned int i;
2161
2162 // Iteratively apply the matrix to the vector
2163 for (i = 1; i <= max_iter; ++i) {
2164
2165 x_prev = x_curr;
2167
2168 // Stop the algorithm when |x_k+1 +- x_k| is
2169 // less then the tolerance in module
2170 if (norm(x_curr - x_prev) <= tolerance ||
2172 break;
2173 }
2174
2175 // The algorithm did not converge
2176 if (i > max_iter) {
2177 TH_MATH_ERROR("eigenpair_power", i, NO_ALGO_CONVERGENCE);
2178 return Type(nan());
2179 }
2180
2181 // Overwrite with the eigenvector
2182 vec_copy(v, x_curr);
2183
2184 return dot(x_curr, A * x_curr);
2185 }
2186
2187
2199 template<typename Matrix, typename Vector>
2201 const Matrix& A, const Vector& x,
2203 unsigned int max_iter = ALGEBRA_EIGEN_ITER) {
2204
2206
2207 if (!is_square(A)) {
2208 TH_MATH_ERROR("eigenvalue_inverse", is_square(A), INVALID_ARGUMENT);
2209 return Type(nan());
2210 }
2211
2212 if (A.rows() != x.size()) {
2213 TH_MATH_ERROR("eigenvalue_inverse", x.size(), INVALID_ARGUMENT);
2214 return Type(nan());
2215 }
2216
2217 // Compute the LU decomposition of A to speed up system solution
2218 Matrix LU = A;
2220
2221 // Compute the first step to initialize the two vectors
2222 Vector x_prev = normalize(x);
2225
2226 // Iteration counter
2227 unsigned int i;
2228
2229 for (i = 1; i <= max_iter; ++i) {
2230
2231 // Solve the linear system using the LU decomposition
2232 x_prev = x_curr;
2234
2235 // Stop the algorithm when |x_k+1 +- x_k| is
2236 // less then the tolerance in module
2237 if (norm(x_curr - x_prev) <= tolerance ||
2239 break;
2240 }
2241
2242 // The algorithm did not converge
2243 if (i > max_iter) {
2244 TH_MATH_ERROR("eigenvalue_inverse", i, NO_ALGO_CONVERGENCE);
2245 return Type(nan());
2246 }
2247
2248 // A and inverse(A) have the same eigenvectors
2249 return dot(x_curr, A * x_curr);
2250 }
2251
2252
2266 template<typename Matrix, typename Vector1, typename Vector2>
2268 const Matrix& A, const Vector1& x, Vector2& v,
2270
2272
2273 if (!is_square(A)) {
2274 TH_MATH_ERROR("eigenpair_inverse", is_square(A), INVALID_ARGUMENT);
2275 return Type(nan());
2276 }
2277
2278 if (A.rows() != x.size()) {
2279 TH_MATH_ERROR("eigenpair_inverse", x.size(), INVALID_ARGUMENT);
2280 return Type(nan());
2281 }
2282
2283 // Compute the LU decomposition of A to speed up system solution
2284 Matrix LU = A;
2286
2287 // Compute the first step to initialize the two vectors
2291
2292 // Iteration counter
2293 unsigned int i;
2294
2295 for (i = 1; i <= max_iter; ++i) {
2296
2297 // Solve the linear system using the LU decomposition
2298 x_prev = x_curr;
2300
2301 // Stop the algorithm when |x_k+1 +- x_k| is
2302 // less then the tolerance in module
2303 if (norm(x_curr - x_prev) <= tolerance ||
2305 break;
2306 }
2307
2308 // The algorithm did not converge
2309 if (i > max_iter) {
2310 TH_MATH_ERROR("eigenpair_inverse", i, NO_ALGO_CONVERGENCE);
2311 return Type(nan());
2312 }
2313
2314 // Overwrite with the eigenvector
2315 vec_copy(v, x_curr);
2316
2317 // A and inverse(A) have the same eigenvectors
2318 return dot(x_curr, A * x_curr);
2319 }
2320
2321
2338 template<typename Matrix, typename Vector, typename T = matrix_element_t<Matrix>>
2340 const Matrix& A, const T& lambda, const Vector& x,
2342
2343 if (!is_square(A)) {
2344 TH_MATH_ERROR("eigenvector_inverse", is_square(A), INVALID_ARGUMENT);
2345 Vector x;
2346 return vec_error(x);
2347 }
2348
2349 if (A.rows() != x.size()) {
2350 TH_MATH_ERROR("eigenvector_inverse", x.size(), INVALID_ARGUMENT);
2351 Vector x;
2352 return vec_error(x);
2353 }
2354
2355 Matrix LU = A;
2357
2358 // Compute the LU decomposition of A
2359 // to speed up system solution
2361
2362 // Compute the first step to initialize the two vectors
2363 Vector v_prev = normalize(x);
2366
2367 // Iteration counter
2368 unsigned int i;
2369
2370 for (i = 1; i <= max_iter; ++i) {
2371
2372 // Solve the linear system using the LU decomposition
2373 v_prev = v_curr;
2375
2376 // Stop the algorithm when |x_k+1 +- x_k| is
2377 // less then the tolerance in module
2378 if (norm(v_curr - v_prev) <= tolerance ||
2380 break;
2381 }
2382
2383 // The algorithm did not converge
2384 if (i > max_iter) {
2385 TH_MATH_ERROR("eigenvector_inverse", i, NO_ALGO_CONVERGENCE);
2386 return vec_error(v_curr);
2387 }
2388
2389 return v_curr;
2390 }
2391
2392
2406 template<typename Matrix, typename Vector, typename T = matrix_element_t<Matrix>>
2408 const Matrix& A, const T& lambda, const Vector& x,
2410
2412
2413 if (!is_square(A)) {
2414 TH_MATH_ERROR("eigenvalue_rayleigh", is_square(A), INVALID_ARGUMENT);
2415 return Type(nan());
2416 }
2417
2418 if (A.rows() != x.size()) {
2419 TH_MATH_ERROR("eigenvalue_rayleigh", x.size(), INVALID_ARGUMENT);
2420 return Type(nan());
2421 }
2422
2423 // Keep track of the shifted matrix
2424 Matrix A_shift = A;
2426
2429
2430 Vector x_prev = normalize(x);
2432
2433 unsigned int i;
2434
2435 for (i = 1; i <= max_iter; ++i) {
2436
2437 // Update the eigenvalue approximation
2438 // using Rayleigh quotient (avoiding normalization)
2441
2442 // Shift the diagonal by the difference between
2443 // subsequent eigenvalues steps, to avoid copying matrix A
2445
2446 // Solve the linear system each time, as the shifted
2447 // matrix is continuously updated
2448 x_prev = x_curr;
2450
2451 // Stop the algorithm when |x_k+1 +- x_k| is
2452 // less then the tolerance in module
2453 if (norm(x_curr - x_prev) <= tolerance ||
2455 break;
2456 }
2457
2458 if (i > max_iter) {
2459 TH_MATH_ERROR("eigenvalue_rayleigh", i, NO_ALGO_CONVERGENCE);
2460 return Type(nan());
2461 }
2462
2463 return dot(x_curr, transform(A, x_curr));
2464 }
2465
2466
2482 template<typename Matrix, typename Vector1, typename Vector2,
2483 typename T = matrix_element_t<Matrix>>
2485 const Matrix& A, const T& lambda, const Vector1& x, Vector2& v,
2487
2489
2490 if (!is_square(A)) {
2491 TH_MATH_ERROR("eigenpair_rayleigh", is_square(A), INVALID_ARGUMENT);
2492 return Type(nan());
2493 }
2494
2495 if (A.rows() != x.size()) {
2496 TH_MATH_ERROR("eigenpair_rayleigh", x.size(), INVALID_ARGUMENT);
2497 return Type(nan());
2498 }
2499
2500 // Keep track of the shifted matrix
2501 Matrix A_shift = A;
2503
2506
2509
2510 unsigned int i;
2511
2512 for (i = 1; i <= max_iter; ++i) {
2513
2514 // Update the eigenvalue approximation
2515 // using Rayleigh quotient (avoiding normalization)
2518
2519 // Shift the diagonal by the difference between
2520 // subsequent eigenvalues steps, to avoid copying matrix A
2522
2523 // Solve the linear system each time, as the shifted
2524 // matrix is continuously updated
2525 x_prev = x_curr;
2527
2528 // Stop the algorithm when |x_k+1 +- x_k| is
2529 // less then the tolerance in module
2530 if (norm(x_curr - x_prev) <= tolerance ||
2532 break;
2533 }
2534
2535 if (i > max_iter) {
2536 TH_MATH_ERROR("eigenpair_rayleigh", i, NO_ALGO_CONVERGENCE);
2537 return Type(nan());
2538 }
2539
2540 // Overwrite with the eigenvector
2541 vec_copy(v, x_curr);
2542
2543 return dot(x_curr, transform(A, x_curr));
2544 }
2545 }
2546}
2547
2548#endif
Complex number in algebraic form .
Definition complex.h:26
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:219
bool is_lower_triangular(const Matrix &m, real tolerance=ALGEBRA_ELEMENT_TOL)
Returns whether the matrix is lower triangular.
Definition algebra.h:1279
Vector1 cross(const Vector1 &v1, const Vector2 &v2)
Compute the cross product between two tridimensional vectors.
Definition algebra.h:373
Matrix2 & mat_sum(Matrix1 &A, const Matrix2 &B)
Sum two matrices and store the result in the first matrix.
Definition algebra.h:764
Matrix3 mat_transpose_mul(const Matrix1 &A, const Matrix2 &B)
Multiply the transpose of a matrix by another matrix, equivalent to the operation .
Definition algebra.h:1049
auto eigenpair_power(const Matrix &A, const Vector1 &x, Vector2 &v, real tolerance=ALGEBRA_EIGEN_TOL, unsigned int max_iter=ALGEBRA_EIGEN_ITER)
Find the biggest eigenvalue in module of a square matrix and its corresponding eigenvector (eigenpai...
Definition algebra.h:2133
Matrix & mat_scalmul(Field a, Matrix &m)
Multiply a matrix by a scalar of any compatible type.
Definition algebra.h:590
bool is_square(const Matrix &m)
Returns whether the matrix is square.
Definition algebra.h:1232
Matrix & invert(Matrix &m)
Invert the given matrix and overwrite it.
Definition algebra.h:1997
auto eigenvalue_power(const Matrix &A, const Vector &x, real tolerance=ALGEBRA_EIGEN_TOL, unsigned int max_iter=ALGEBRA_EIGEN_ITER)
Find the biggest eigenvalue in module of a square matrix using the power method (Von Mises iteration...
Definition algebra.h:2072
auto trace(const Matrix &m)
Compute the trace of the given matrix.
Definition algebra.h:555
auto eigenpair_inverse(const Matrix &A, const Vector1 &x, Vector2 &v, real tolerance=ALGEBRA_EIGEN_TOL, unsigned int max_iter=ALGEBRA_EIGEN_ITER)
Find the eigenvalue with the smallest inverse of a square matrix and its corresponding eigenvector,...
Definition algebra.h:2267
Matrix & make_identity(Matrix &m)
Overwrite a matrix with the identity matrix.
Definition algebra.h:77
Matrix2 & mat_diff(Matrix1 &A, const Matrix2 &B)
Subtract two matrices and store the result in the first matrix.
Definition algebra.h:832
Matrix & mat_swap_cols(Matrix &A, unsigned int col1, unsigned int col2)
Swap two columns of a matrix, given the matrix and the two indices of the columns.
Definition algebra.h:197
Vector & make_normalized(Vector &v)
Normalize a given vector overwriting it.
Definition algebra.h:327
Vector1 & vec_copy(Vector1 &dest, const Vector2 &src)
Copy a vector by overwriting another.
Definition algebra.h:143
Vector & solve_lu_inplace(const Matrix &A, Vector &b)
Solve the linear system , finding , where the matrix A has already undergone in-place LU decompositio...
Definition algebra.h:1681
bool is_symmetric(const Matrix &m, real tolerance=ALGEBRA_ELEMENT_TOL)
Returns whether the matrix is symmetric.
Definition algebra.h:1259
Vector solve(const Matrix &A, const Vector &b)
Solve the linear system , finding using the best available algorithm.
Definition algebra.h:1867
bool mat_equals(const Matrix1 &A, const Matrix2 &B, real tolerance=10 *MACH_EPSILON)
Checks whether two matrices are equal.
Definition algebra.h:1109
Vector transform(const Matrix &A, const Vector &v)
Returns the matrix transformation of a vector.
Definition algebra.h:702
bool is_upper_triangular(const Matrix &m, real tolerance=ALGEBRA_ELEMENT_TOL)
Returns whether the matrix is upper triangular.
Definition algebra.h:1299
bool is_diagonal(const Matrix &m, real tolerance=10 *MACH_EPSILON)
Returns whether the matrix is diagonal.
Definition algebra.h:1243
Vector & apply_transform(const Matrix &A, Vector &v)
Apply a matrix transformation to a vector and store the result in the vector.
Definition algebra.h:675
real density(const Matrix &A, real tolerance=1E-12)
Compute the density of a matrix, counting the proportion of non-zero (bigger in module than the given...
Definition algebra.h:1323
Vector eigenvector_inverse(const Matrix &A, const T &lambda, const Vector &x, real tolerance=ALGEBRA_EIGEN_TOL, unsigned int max_iter=ALGEBRA_EIGEN_ITER)
Find the eigenvector associated with a given approximated eigenvalue using the inverse power method.
Definition algebra.h:2339
void decompose_lu(const Matrix1 &A, Matrix2 &L, Matrix3 &U)
Decompose a square matrix to two triangular matrices, L and U where L is lower and U is upper,...
Definition algebra.h:1363
Matrix & mat_swap_rows(Matrix &A, unsigned int row1, unsigned int row2)
Swap two rows of a matrix, given the matrix and the two indices of the rows.
Definition algebra.h:161
Vector solve_lu(Matrix A, Vector b)
Solve the linear system , finding .
Definition algebra.h:1729
real sparsity(const Matrix &A, real tolerance=1E-12)
Compute the sparsity of a matrix, counting the proportion of zero (smaller in module than the given t...
Definition algebra.h:1346
Vector solve_triangular(const Matrix &T, const Vector &b)
Solve the linear system for triangular .
Definition algebra.h:1655
Matrix & mat_error(Matrix &m)
Overwrite the given matrix with the error matrix with NaN values on the diagonal and zeroes everywher...
Definition algebra.h:44
Matrix decompose_cholesky(const Matrix &A)
Decompose a symmetric positive definite matrix into a triangular matrix so that using Cholesky decom...
Definition algebra.h:1462
auto dot(const Vector1 &v, const Vector2 &w)
Computes the dot product between two vectors, using the Hermitian form if needed.
Definition algebra.h:351
Matrix3 mat_mul(const Matrix1 &A, const Matrix2 &B)
Multiply two matrices and store the result in the first matrix, equivalent to the operation .
Definition algebra.h:977
Matrix & mat_zeroes(Matrix &m)
Overwrite a matrix with all zeroes.
Definition algebra.h:93
auto pair_inner_product(const Type &v_i, const Type &w_i)
Compute the contribution of the inner product between a pair of elements of two vectors,...
Definition algebra.h:252
auto diagonal_product(const Matrix &m)
Compute the product of the elements of the main diagonal of a generic matrix.
Definition algebra.h:573
auto eigenvalue_rayleigh(const Matrix &A, const T &lambda, const Vector &x, real tolerance=ALGEBRA_EIGEN_TOL, unsigned int max_iter=ALGEBRA_EIGEN_ITER)
Compute an eigenvalue of a square matrix using the Rayleigh quotient iteration method.
Definition algebra.h:2407
auto eigenpair_rayleigh(const Matrix &A, const T &lambda, const Vector1 &x, Vector2 &v, real tolerance=ALGEBRA_EIGEN_TOL, unsigned int max_iter=ALGEBRA_EIGEN_ITER)
Compute an eigenvalue of a square matrix and its corresponding eigenvector (eigenpair) using the Rayl...
Definition algebra.h:2484
auto rayleigh_quotient(const Matrix &A, const Vector &x)
Compute the Rayleigh quotient of a vector with respect to a matrix.
Definition algebra.h:2046
auto sqr_norm(const Vector &v)
Returns the square of the Euclidean/Hermitian norm of the given vector.
Definition algebra.h:275
auto det(const Matrix &A)
Compute the determinant of a square matrix.
Definition algebra.h:2024
Matrix & mat_shift_diagonal(Matrix &A, const Type &sigma)
Shift the diagonal of a matrix by the given amount, overwriting the matrix itself,...
Definition algebra.h:233
Vector solve_triangular_lower(const Matrix &L, const Vector &b)
Solve the linear system for lower triangular .
Definition algebra.h:1564
Matrix2 & mat_lincomb(Field1 alpha, Matrix1 &A, Field2 beta, const Matrix2 &B)
Compute the linear combination of two matrices and store the result in the first matrix.
Definition algebra.h:901
Matrix & decompose_cholesky_inplace(Matrix &A)
Decompose a symmetric positive definite matrix in-place, overwriting the starting matrix,...
Definition algebra.h:1518
Vector2 & vec_sum(Vector1 &v1, const Vector2 &v2)
Sum two vectors and store the result in the first vector.
Definition algebra.h:1135
Matrix3 mat_mul_transpose(const Matrix1 &A, const Matrix2 &B)
Multiply a matrix by the transpose of another matrix, equivalent to the operation .
Definition algebra.h:1081
Matrix1 & mat_copy(Matrix1 &dest, const Matrix2 &src)
Copy a matrix by overwriting another.
Definition algebra.h:125
Vector2 & vec_diff(Vector1 &v1, const Vector2 &v2)
Subtract two vectors and store the result in the first vector.
Definition algebra.h:1183
MatrixT hermitian(const Matrix &m)
Returns the hermitian of the given matrix.
Definition algebra.h:478
Matrix & make_transposed(Matrix &m)
Transpose the given matrix.
Definition algebra.h:420
Vector solve_triangular_upper(const Matrix &U, const Vector &b)
Solve the linear system for upper triangular .
Definition algebra.h:1606
Vector & vec_zeroes(Vector &v)
Overwrite a vector with all zeroes.
Definition algebra.h:109
Vector & vec_error(Vector &v)
Overwrite the given vector with the error vector with NaN values.
Definition algebra.h:62
Vector solve_cholesky(const Matrix &L, const Vector &b)
Solve a linear system defined by a symmetric positive definite matrix, using the Cholesky decomposit...
Definition algebra.h:1818
Matrix1 & inverse(Matrix1 &dest, const Matrix2 &src)
Invert the given matrix.
Definition algebra.h:1883
auto norm(const Vector &v)
Returns the Euclidean/Hermitian norm of the given vector.
Definition algebra.h:292
Vector normalize(const Vector &v)
Returns the normalized vector.
Definition algebra.h:302
auto eigenvalue_inverse(const Matrix &A, const Vector &x, real tolerance=ALGEBRA_EIGEN_TOL, unsigned int max_iter=ALGEBRA_EIGEN_ITER)
Find the eigenvalue with the smallest inverse of a square matrix, using the inverse power method wit...
Definition algebra.h:2200
Matrix & decompose_lu_inplace(Matrix &A)
Decompose a square matrix to two triangular matrices, L and U where L is lower and U is upper,...
Definition algebra.h:1432
Matrix & make_hermitian(Matrix &m)
Compute the hermitian of a given matrix and overwrite it.
Definition algebra.h:497
Vector & vec_scalmul(Field a, Vector &v)
Multiply a vector by a scalar of any compatible type.
Definition algebra.h:634
MatrixT transpose(const Matrix &m)
Returns the transpose of the given matrix.
Definition algebra.h:403
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
double real
A real number, defined as a floating point type.
Definition constants.h:198
auto min(const Vector &X)
Finds the minimum value inside a dataset.
Definition dataset.h:351
dual2 sqrt(dual2 x)
Compute the square root of a second order dual number.
Definition dual2_functions.h:54
bool is_nan(const T &x)
Check whether a generic variable is (equivalent to) a NaN number.
Definition error.h:70
dual2 abs(dual2 x)
Compute the absolute value of a second order dual number.
Definition dual2_functions.h:198
std::remove_reference_t< decltype(std::declval< Structure >()(0, 0))> matrix_element_t
Extract the type of a matrix (or any doubly indexable container) from its operator().
Definition core_traits.h:140
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
constexpr real ALGEBRA_ELEMENT_TOL
Tolerance for the elements of matrices.
Definition constants.h:270
auto sum(const Vector &X)
Compute the sum of a vector of real values using pairwise summation to reduce round-off error.
Definition dataset.h:219
dual2 conjugate(dual2 x)
Return the conjugate of a second order dual number.
Definition dual2_functions.h:35
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition constants.h:207
constexpr real E
The Euler mathematical constant (e)
Definition constants.h:237
constexpr real ALGEBRA_EIGEN_ITER
Maximum number of iterations for eigensolvers.
Definition constants.h:276
real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
constexpr real ALGEBRA_EIGEN_TOL
Tolerance for eigensolvers.
Definition constants.h:273