Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
Matrix.cpp
Go to the documentation of this file.
1// #include "Matrix.h"
2
3// std::ostream &operator<<(std::ostream &os, RealMatrix A) {
4// os << "[";
5// for (size_t i = 0; i < A.rows(); i++) {
6// os << (i == 0 ? "[ " : " [ ");
7// for (size_t j = 0; j < A.cols(); j++) {
8// os << A.at(i, j) << (j == A.cols() - 1 ? " " : ", ");
9// }
10// if (i != A.rows() - 1)
11// os << "]" << std::endl;
12// }
13// os << "]]" << std::endl;
14// return os;
15// }
16
17// RealMatrix eye(std::size_t n)
18// {
19// std::vector<double> new_data(n * n, 0.0);
20
21// for (size_t i = 0; i < n; ++i) {
22// new_data[i * (n + 1)] = 1.0;
23// }
24
25// return RealMatrix(new_data, n, n);
26// }
27
28// RealMatrix diag(const gsl_vector *X) {
29// const std::size_t n = X->size;
30// std::vector<double> new_data(n * n, 0.0);
31
32// for (size_t i = 0; i < n; ++i) {
33// new_data[i * (n + 1)] = gsl_vector_get(X, i);
34// }
35
36// return RealMatrix(new_data, n, n);
37// }
38
39// RealMatrix nearest_psd(RealMatrix R, double thr) {
40// if (R.rows() != R.cols())
41// throw std::invalid_argument("Matrix should be square");
42
43// size_t n = R.rows();
44
45// // Symmetrize
46// R = 0.5 * (R + R.transpose());
47
48// // Compute eigensystem, clip eigenvalues and recombine
49// EigenSystem e = R.eig();
50// for (size_t i = 0; i < n; i++) {
51// e.D.at(i, i) = std::max(e.D.at(i, i), thr);
52// }
53// RealMatrix R_psd = e.P * e.D * e.P.transpose();
54
55// // Normalize, re-symmetrize and enforce unit correlations on the diagonal
56// std::vector<double> inv_sqrt_diag(n);
57// for (size_t i = 0; i < n; ++i)
58// inv_sqrt_diag[i] = 1.0 / std::sqrt(R_psd.at(i, i));
59
60// for (size_t i = 0; i < n; ++i)
61// for (size_t j = 0; j < n; ++j)
62// R_psd.at(i, j) *= inv_sqrt_diag[i] * inv_sqrt_diag[j];
63
64// R_psd = 0.5 * (R_psd + R_psd.transpose());
65
66// for (size_t i = 0; i < n; i++)
67// R_psd.at(i, i) = 1.0;
68
69// return R_psd;
70// }
71
72// RealMatrix cholesky_L(RealMatrix R) {
73// if (R.rows() != R.cols())
74// throw std::invalid_argument("Matrix should be square");
75
76// gsl_matrix_sptr R_gsl = R.to_gsl_matrix();
77
78// if (gsl_linalg_cholesky_decomp1(R_gsl.get()) != GSL_SUCCESS)
79// throw std::runtime_error("Cholesky decomposition failed (matrix not PSD)");
80
81// const size_t n = R.rows();
82// RealMatrix L(n, n);
83
84// for (size_t i = 0; i < n; ++i) {
85// for (size_t j = 0; j <= i; ++j) {
86// L.at(i, j) = gsl_matrix_get(R_gsl.get(), i, j);
87// }
88// }
89
90// return L;
91// }
92
93// RealMatrix block_diag(const std::vector<RealMatrix> &blocks) {
94// if (blocks.empty())
95// return RealMatrix{};
96
97// std::size_t total_dim = 0;
98
99// for (const auto& B : blocks) {
100// if (B.rows() != B.cols())
101// throw std::invalid_argument("block_diag: matrices must be square");
102// total_dim += B.rows();
103// }
104
105// RealMatrix M(total_dim, total_dim);
106
107// std::size_t offset = 0;
108// for (const auto& B : blocks) {
109// const std::size_t n = B.rows();
110// for (std::size_t i = 0; i < n; ++i)
111// for (std::size_t j = 0; j < n; ++j)
112// M.unchecked_at(offset + i, offset + j) = B.unchecked_at(i, j);
113// offset += n;
114// }
115
116// return M;
117// }
118
119// RealMatrix::RealMatrix(std::vector<double> data_, std::size_t rows, std::size_t cols) : data(data_), rows_(rows), cols_(cols) {
120// if (data.size() != rows * cols)
121// throw std::invalid_argument("Data size does not match matrix shape.");
122// }
123
124// RealMatrix::RealMatrix(std::vector<std::vector<double>> data) {
125// if (data.size() == 0)
126// throw std::invalid_argument("Matrix should not be sizeless.");
127
128// this->rows_ = data.size();
129
130// std::size_t size_2 = data[0].size();
131// for (auto& row : data) {
132// if (row.size() != size_2)
133// throw std::invalid_argument("All rows should have the same number of columns.");
134// }
135
136// this->cols_ = size_2;
137
138// std::vector<double> new_data (rows_ * cols_, 0.0);
139// for (size_t i = 0; i < rows_; i++) {
140// for (size_t j = 0; j < cols_; j++) {
141// new_data[i * cols_ + j] = data[i][j];
142// }
143// }
144
145// this->data = std::move(new_data);
146// }
147
148// RealMatrix::RealMatrix(std::size_t rows, std::size_t cols) : data(rows * cols, 0.0), rows_(rows), cols_(cols)
149// {}
150
151// double &RealMatrix::at(size_t i, size_t j) {
152// if (i >= rows_ || j >= cols_)
153// throw std::out_of_range("Indices outside matrix shape");
154
155// return this->data[i * cols_ + j];
156// }
157
158// const double &RealMatrix::at(size_t i, size_t j) const {
159// if (i >= rows_ || j >= cols_)
160// throw std::out_of_range("Indices outside matrix shape");
161
162// return this->data[i * cols_ + j];
163// }
164
165// std::size_t RealMatrix::rows() const {
166// return this->rows_;
167// }
168
169// std::size_t RealMatrix::cols() const {
170// return this->cols_;
171// }
172
173// void RealMatrix::remove_row(std::size_t row_idx) {
174// std::vector<double> new_data = std::vector<double>((this->rows_ - 1) * this->cols_, 0.0);
175
176// for (size_t i = 0; i < this->rows_; i++) {
177// if (i == row_idx) continue;
178// std::size_t new_i = i < row_idx ? i : i - 1;
179// for (size_t j = 0; j < this->cols_; j++) {
180// new_data[new_i * this->cols_ + j] = this->at(i, j);
181// }
182// }
183
184// this->data = std::move(new_data);
185// this->rows_--;
186// }
187
188// void RealMatrix::remove_column(std::size_t col_idx) {
189// std::vector<double> new_data = std::vector<double>(this->rows_ * (this->cols_ - 1), 0.0);
190
191// for (size_t j = 0; j < this->cols_; j++) {
192// if (j == col_idx) continue;
193// std::size_t new_j = j < col_idx ? j : j - 1;
194// for (size_t i = 0; i < this->rows_; i++) {
195// new_data[i * (this->cols_ - 1) + new_j] = this->at(i, j);
196// }
197// }
198
199// this->data = std::move(new_data);
200// this->cols_--;
201// }
202
203// void RealMatrix::remove_row_and_column(std::size_t dim_idx) {
204// std::vector<double> new_data((this->rows_ - 1) * (this->cols_ - 1), 0.0);
205
206// for (size_t i = 0; i < this->rows_; i++) {
207// if (i == dim_idx) continue;
208// std::size_t new_i = i < dim_idx ? i : i - 1;
209
210// for (size_t j = 0; j < this->cols_; j++) {
211// if (j == dim_idx) continue;
212// std::size_t new_j = j < dim_idx ? j : j - 1;
213// new_data[new_i * (this->cols_ - 1) + new_j] = this->at(i, j);
214// }
215// }
216
217// this->data = std::move(new_data);
218// this->rows_--;
219// this->cols_--;
220// }
221
222// RealMatrix RealMatrix::from_gsl_copy(const gsl_matrix* A) {
223// const size_t n = A->size1;
224// const size_t m = A->size2;
225// std::vector<double> data(n * m, 0.0);
226
227// for (size_t i = 0; i < n; ++i) {
228// for (size_t j = 0; j < m; ++j) {
229// data[i * m + j] = gsl_matrix_get(A, i, j);
230// }
231// }
232
233// return RealMatrix(std::move(data), n, m);
234// }
235
236// gsl_matrix_sptr RealMatrix::to_gsl_matrix() const {
237// gsl_matrix_sptr A = make_gsl_matrix(rows_, cols_);
238
239// for (size_t i = 0; i < rows_; ++i) {
240// for (size_t j = 0; j < cols_; ++j) {
241// gsl_matrix_set(A.get(), i, j, this->data[i * cols_ + j]);
242// }
243// }
244
245// return A;
246// }
247
248// bool RealMatrix::is_symmetric() const {
249// if (rows_ != cols_)
250// return false;
251
252// for (std::size_t i = 0; i < rows_; i++) {
253// for (std::size_t j = i + 1; j < cols_; j++) {
254// if (!fpeq(data[i * cols_ + j], data[j * cols_ + i]))
255// return false;
256// }
257// }
258
259// return true;
260// }
261
262// RealMatrix RealMatrix::transpose() const {
263// auto new_data = std::vector(this->cols_ * this->rows_, 0.0);
264
265// for (size_t i = 0; i < this->rows_; ++i) {
266// for (size_t j = 0; j < this->cols_; ++j) {
267// new_data[j * rows_ + i] = this->unchecked_at(i, j);
268// }
269// }
270
271// return RealMatrix(std::move(new_data), cols_, rows_);
272// }
273
274// EigenSystem RealMatrix::eig() const {
275// if (rows_ != cols_)
276// throw std::invalid_argument("Eigen decomposition requires a square matrix.");
277
278// if (!is_symmetric())
279// throw std::runtime_error("Matrix is not symmetric.");
280
281// gsl_matrix_sptr M = this->to_gsl_matrix();
282// gsl_vector_sptr eval = make_gsl_vector (rows_);
283// gsl_matrix_sptr evec = make_gsl_matrix (rows_, rows_);
284// gsl_eigen_workspace_sptr w = make_eigen_workspace(rows_);
285// gsl_eigen_symmv (M.get(), eval.get(), evec.get(), w.get());
286// gsl_eigen_symmv_sort(eval.get(), evec.get(), GSL_EIGEN_SORT_VAL_DESC);
287
288// EigenSystem e;
289// e.P = from_gsl_copy(evec.get());
290// e.D = diag(eval.get());
291
292// return e;
293// }
294
295// SignedLogDet RealMatrix::slogdet() const {
296// if (rows_ != cols_)
297// throw std::invalid_argument("LU decomposition requires a square matrix.");
298
299// gsl_matrix_sptr M = this->to_gsl_matrix();
300// gsl_permutation_sptr p = make_gsl_permutation(rows_);
301// int p_sign;
302
303// if (gsl_linalg_LU_decomp(M.get(), p.get(), &p_sign) != GSL_SUCCESS) {
304// return {0.0, 1};
305// }
306
307// SignedLogDet sld;
308// sld.logdet = gsl_linalg_LU_lndet(M.get());
309// sld.sign = gsl_linalg_LU_sgndet(M.get(), p_sign);
310
311// return sld;
312// }
313
314// RealMatrix RealMatrix::inv() const {
315// if (rows_ != cols_)
316// throw std::invalid_argument("inversion requires a square matrix.");
317
318// gsl_matrix_sptr M = this->to_gsl_matrix();
319// gsl_matrix_sptr M_inv = make_gsl_matrix(rows_, cols_);
320// gsl_permutation_sptr p = make_gsl_permutation(rows_);
321// int p_sign;
322
323// if (gsl_linalg_LU_decomp(M.get(), p.get(), &p_sign) != GSL_SUCCESS) {
324// throw std::runtime_error("Matrix is singular");
325// }
326
327// gsl_linalg_LU_invert(M.get(), p.get(), M_inv.get());
328
329// return from_gsl_copy(M_inv.get());
330// }
331
332// RealMatrix RealMatrix::operator-() const {
333// auto new_data = std::vector(this->rows_ * this->cols_, 0.0);
334
335// for (size_t i = 0; i < this->rows_; ++i) {
336// for (size_t j = 0; j < this->cols_; ++j) {
337// new_data[i * cols_ + j] = -this->unchecked_at(i, j);
338// }
339// }
340
341// return RealMatrix(std::move(new_data), rows_, cols_);
342// }
343
344// RealMatrix &RealMatrix::operator+=(const RealMatrix &rhs) {
345// if (this->rows_ != rhs.rows() || this->cols_ != rhs.cols())
346// throw std::invalid_argument("Matrices should have the same shape.");
347
348// for (size_t i = 0; i < this->rows_; ++i) {
349// for (size_t j = 0; j < this->cols_; ++j) {
350// this->unchecked_at(i, j) += rhs.unchecked_at(i, j);
351// }
352// }
353
354// return *this;
355// }
356
357// RealMatrix &RealMatrix::operator-=(const RealMatrix &rhs) {
358// if (this->rows_ != rhs.rows() || this->cols_ != rhs.cols())
359// throw std::invalid_argument("Matrices should have the same shape.");
360
361// for (size_t i = 0; i < this->rows_; ++i) {
362// for (size_t j = 0; j < this->cols_; ++j) {
363// this->unchecked_at(i, j) -= rhs.unchecked_at(i, j);
364// }
365// }
366
367// return *this;
368// }
369
370// RealMatrix &RealMatrix::operator*=(const RealMatrix &rhs) {
371// if (this->cols_ != rhs.rows())
372// throw std::invalid_argument("Matrices don't have the right shape to be multiplied.");
373
374// const size_t M = rows_;
375// const size_t N = rhs.cols_;
376// const size_t K = cols_;
377
378// std::vector<double> result(M * N, 0.0);
379
380// constexpr size_t BLOCK = 32;
381
382// for (size_t ii = 0; ii < M; ii += BLOCK) {
383// for (size_t kk = 0; kk < K; kk += BLOCK) {
384// for (size_t jj = 0; jj < N; jj += BLOCK) {
385
386// size_t i_max = std::min(ii + BLOCK, M);
387// size_t k_max = std::min(kk + BLOCK, K);
388// size_t j_max = std::min(jj + BLOCK, N);
389
390// for (size_t i = ii; i < i_max; ++i) {
391// for (size_t k = kk; k < k_max; ++k) {
392// double aik = unchecked_at(i, k);
393// const double* rhs_row = &rhs.data[k * rhs.cols_];
394
395// double* res_row = &result[i * N];
396
397// for (size_t j = jj; j < j_max; ++j) {
398// res_row[j] += aik * rhs_row[j];
399// }
400// }
401// }
402// }
403// }
404// }
405
406// data = std::move(result);
407// cols_ = N;
408// return *this;
409// }
410
411// RealMatrix &RealMatrix::operator*=(double scalar) {
412// for (double& x: data)
413// x *= scalar;
414
415// return *this;
416// }
417
418// RealMatrix &RealMatrix::operator/=(double scalar) {
419// if (fpeq(scalar, 0.0))
420// throw std::invalid_argument("Division by zero.");
421
422// for (double& x: data)
423// x /= scalar;
424
425// return *this;
426// }
427
428// double &RealMatrix::unchecked_at(size_t i, size_t j) {
429// return data[i * cols_ + j];
430// }
431
432// const double &RealMatrix::unchecked_at(size_t i, size_t j) const {
433// return data[i * cols_ + j];
434// }
435
436#include "Matrix.h"
437#include <gsl/gsl_errno.h>
438
439std::ostream &operator<<(std::ostream &os, RealMatrix A) {
440 os << "[";
441 for (size_t i = 0; i < A.rows(); i++) {
442 os << (i == 0 ? "[ " : " [ ");
443 for (size_t j = 0; j < A.cols(); j++) {
444 os << A.at(i, j) << (j == A.cols() - 1 ? " " : ", ");
445 }
446 if (i != A.rows() - 1)
447 os << "]" << std::endl;
448 }
449 os << "]]" << std::endl;
450 return os;
451}
452
453RealMatrix eye(std::size_t n)
454{
455 std::vector<double> new_data(n * n, 0.0);
456
457 for (size_t i = 0; i < n; ++i) {
458 new_data[i * (n + 1)] = 1.0;
459 }
460
461 return RealMatrix(new_data, n, n);
462}
463
464RealMatrix diag(const gsl_vector *X) {
465 const std::size_t n = X->size;
466 std::vector<double> new_data(n * n, 0.0);
467
468 for (size_t i = 0; i < n; ++i) {
469 new_data[i * (n + 1)] = gsl_vector_get(X, i);
470 }
471
472 return RealMatrix(new_data, n, n);
473}
474
475RealMatrix nearest_psd(RealMatrix R, double thr) {
476 if (R.rows() != R.cols())
477 throw std::invalid_argument("Matrix should be square");
478
479 size_t n = R.rows();
480
481 // Symmetrize
482 R = 0.5 * (R + R.transpose());
483
484 // Compute eigensystem, clip eigenvalues and recombine
485 EigenSystem e = R.eig();
486 for (size_t i = 0; i < n; i++) {
487 e.D.at(i, i) = std::max(e.D.at(i, i), thr);
488 }
489 RealMatrix R_psd = e.P * e.D * e.P.transpose();
490
491 // Normalize, re-symmetrize and enforce unit correlations on the diagonal
492 std::vector<double> inv_sqrt_diag(n);
493 for (size_t i = 0; i < n; ++i)
494 inv_sqrt_diag[i] = 1.0 / std::sqrt(R_psd.at(i, i));
495
496 for (size_t i = 0; i < n; ++i)
497 for (size_t j = 0; j < n; ++j)
498 R_psd.at(i, j) *= inv_sqrt_diag[i] * inv_sqrt_diag[j];
499
500 R_psd = 0.5 * (R_psd + R_psd.transpose());
501
502 for (size_t i = 0; i < n; i++)
503 R_psd.at(i, i) = 1.0;
504
505 return R_psd;
506}
507
508RealMatrix cholesky_L(RealMatrix R) {
509 if (R.rows() != R.cols())
510 throw std::invalid_argument("Matrix should be square");
511
512 gsl_matrix_sptr R_gsl = R.to_gsl_matrix();
513
514 if (gsl_linalg_cholesky_decomp1(R_gsl.get()) != GSL_SUCCESS)
515 throw std::runtime_error("Cholesky decomposition failed (matrix not PSD)");
516
517 const size_t n = R.rows();
518 RealMatrix L(n, n);
519
520 for (size_t i = 0; i < n; ++i) {
521 for (size_t j = 0; j <= i; ++j) {
522 L.at(i, j) = gsl_matrix_get(R_gsl.get(), i, j);
523 }
524 }
525
526 return L;
527}
528
529RealMatrix block_diag(const std::vector<RealMatrix> &blocks) {
530 if (blocks.empty())
531 return RealMatrix{};
532
533 std::size_t total_dim = 0;
534
535 for (const auto& B : blocks) {
536 if (B.rows() != B.cols())
537 throw std::invalid_argument("block_diag: matrices must be square");
538 total_dim += B.rows();
539 }
540
541 RealMatrix M(total_dim, total_dim);
542
543 std::size_t offset = 0;
544 for (const auto& B : blocks) {
545 const std::size_t n = B.rows();
546 for (std::size_t i = 0; i < n; ++i)
547 for (std::size_t j = 0; j < n; ++j)
548 M.unchecked_at(offset + i, offset + j) = B.unchecked_at(i, j);
549 offset += n;
550 }
551
552 return M;
553}
554
555RealMatrix::RealMatrix(std::vector<double> data_, std::size_t rows, std::size_t cols) : data(data_), rows_(rows), cols_(cols) {
556 if (data.size() != rows * cols)
557 throw std::invalid_argument("Data size does not match matrix shape.");
558}
559
560RealMatrix::RealMatrix(std::vector<std::vector<double>> data) {
561 if (data.size() == 0)
562 throw std::invalid_argument("Matrix should not be sizeless.");
563
564 this->rows_ = data.size();
565
566 std::size_t size_2 = data[0].size();
567 for (auto& row : data) {
568 if (row.size() != size_2)
569 throw std::invalid_argument("All rows should have the same number of columns.");
570 }
571
572 this->cols_ = size_2;
573
574 std::vector<double> new_data (rows_ * cols_, 0.0);
575 for (size_t i = 0; i < rows_; i++) {
576 for (size_t j = 0; j < cols_; j++) {
577 new_data[i * cols_ + j] = data[i][j];
578 }
579 }
580
581 this->data = std::move(new_data);
582}
583
584RealMatrix::RealMatrix(std::size_t rows, std::size_t cols) : data(rows * cols, 0.0), rows_(rows), cols_(cols)
585{}
586
587double &RealMatrix::at(size_t i, size_t j) {
588 if (i >= rows_ || j >= cols_)
589 throw std::out_of_range("Indices outside matrix shape");
590
591 return this->data[i * cols_ + j];
592}
593
594const double &RealMatrix::at(size_t i, size_t j) const {
595 if (i >= rows_ || j >= cols_)
596 throw std::out_of_range("Indices outside matrix shape");
597
598 return this->data[i * cols_ + j];
599}
600
601std::size_t RealMatrix::rows() const {
602 return this->rows_;
603}
604
605std::size_t RealMatrix::cols() const {
606 return this->cols_;
607}
608
609void RealMatrix::remove_row(std::size_t row_idx) {
610 std::vector<double> new_data = std::vector<double>((this->rows_ - 1) * this->cols_, 0.0);
611
612 for (size_t i = 0; i < this->rows_; i++) {
613 if (i == row_idx) continue;
614 std::size_t new_i = i < row_idx ? i : i - 1;
615 for (size_t j = 0; j < this->cols_; j++) {
616 new_data[new_i * this->cols_ + j] = this->at(i, j);
617 }
618 }
619
620 this->data = std::move(new_data);
621 this->rows_--;
622}
623
624void RealMatrix::remove_column(std::size_t col_idx) {
625 std::vector<double> new_data = std::vector<double>(this->rows_ * (this->cols_ - 1), 0.0);
626
627 for (size_t j = 0; j < this->cols_; j++) {
628 if (j == col_idx) continue;
629 std::size_t new_j = j < col_idx ? j : j - 1;
630 for (size_t i = 0; i < this->rows_; i++) {
631 new_data[i * (this->cols_ - 1) + new_j] = this->at(i, j);
632 }
633 }
634
635 this->data = std::move(new_data);
636 this->cols_--;
637}
638
639void RealMatrix::remove_row_and_column(std::size_t dim_idx) {
640 std::vector<double> new_data((this->rows_ - 1) * (this->cols_ - 1), 0.0);
641
642 for (size_t i = 0; i < this->rows_; i++) {
643 if (i == dim_idx) continue;
644 std::size_t new_i = i < dim_idx ? i : i - 1;
645
646 for (size_t j = 0; j < this->cols_; j++) {
647 if (j == dim_idx) continue;
648 std::size_t new_j = j < dim_idx ? j : j - 1;
649 new_data[new_i * (this->cols_ - 1) + new_j] = this->at(i, j);
650 }
651 }
652
653 this->data = std::move(new_data);
654 this->rows_--;
655 this->cols_--;
656}
657
658RealMatrix RealMatrix::from_gsl_copy(const gsl_matrix* A) {
659 const size_t n = A->size1;
660 const size_t m = A->size2;
661 std::vector<double> data(n * m, 0.0);
662
663 for (size_t i = 0; i < n; ++i) {
664 for (size_t j = 0; j < m; ++j) {
665 data[i * m + j] = gsl_matrix_get(A, i, j);
666 }
667 }
668
669 return RealMatrix(std::move(data), n, m);
670}
671
672gsl_matrix_sptr RealMatrix::to_gsl_matrix() const {
673 gsl_matrix_sptr A = make_gsl_matrix(rows_, cols_);
674
675 for (size_t i = 0; i < rows_; ++i) {
676 for (size_t j = 0; j < cols_; ++j) {
677 gsl_matrix_set(A.get(), i, j, this->data[i * cols_ + j]);
678 }
679 }
680
681 return A;
682}
683
684bool RealMatrix::is_symmetric() const {
685 if (rows_ != cols_)
686 return false;
687
688 for (std::size_t i = 0; i < rows_; i++) {
689 for (std::size_t j = i + 1; j < cols_; j++) {
690 if (!fpeq(data[i * cols_ + j], data[j * cols_ + i]))
691 return false;
692 }
693 }
694
695 return true;
696}
697
698RealMatrix RealMatrix::transpose() const {
699 auto new_data = std::vector(this->cols_ * this->rows_, 0.0);
700
701 for (size_t i = 0; i < this->rows_; ++i) {
702 for (size_t j = 0; j < this->cols_; ++j) {
703 new_data[j * rows_ + i] = this->unchecked_at(i, j);
704 }
705 }
706
707 return RealMatrix(std::move(new_data), cols_, rows_);
708}
709
710EigenSystem RealMatrix::eig() const {
711 if (rows_ != cols_)
712 throw std::invalid_argument("Eigen decomposition requires a square matrix.");
713
714 if (!is_symmetric())
715 throw std::runtime_error("Matrix is not symmetric.");
716
717 gsl_matrix_sptr M = this->to_gsl_matrix();
718 gsl_vector_sptr eval = make_gsl_vector (rows_);
719 gsl_matrix_sptr evec = make_gsl_matrix (rows_, rows_);
720 gsl_eigen_workspace_sptr w = make_eigen_workspace(rows_);
721 gsl_eigen_symmv (M.get(), eval.get(), evec.get(), w.get());
722 gsl_eigen_symmv_sort(eval.get(), evec.get(), GSL_EIGEN_SORT_VAL_DESC);
723
724 EigenSystem e;
725 e.P = from_gsl_copy(evec.get());
726 e.D = diag(eval.get());
727
728 return e;
729}
730
731SignedLogDet RealMatrix::slogdet() const {
732 if (rows_ != cols_)
733 throw std::invalid_argument("LU decomposition requires a square matrix.");
734
735 gsl_matrix_sptr M = this->to_gsl_matrix();
736 gsl_permutation_sptr p = make_gsl_permutation(rows_);
737 int p_sign;
738
739 if (gsl_linalg_LU_decomp(M.get(), p.get(), &p_sign) != GSL_SUCCESS) {
740 return {0.0, 1};
741 }
742
743 SignedLogDet sld;
744 sld.logdet = gsl_linalg_LU_lndet(M.get());
745 sld.sign = gsl_linalg_LU_sgndet(M.get(), p_sign);
746
747 return sld;
748}
749
750RealMatrix RealMatrix::inv() const {
751 if (rows_ != cols_)
752 throw std::invalid_argument("inversion requires a square matrix.");
753
754 gsl_error_handler_t* old_handler = gsl_set_error_handler_off();
755
756 gsl_matrix_sptr M = this->to_gsl_matrix();
757 gsl_matrix_sptr M_inv = make_gsl_matrix(rows_, cols_);
758 gsl_permutation_sptr p = make_gsl_permutation(rows_);
759 int p_sign = 0;
760
761 const int decomp_status = gsl_linalg_LU_decomp(M.get(), p.get(), &p_sign);
762 if (decomp_status != GSL_SUCCESS) {
763 gsl_set_error_handler(old_handler);
764 throw std::runtime_error("Matrix LU decomposition failed");
765 }
766
767 const int invert_status = gsl_linalg_LU_invert(M.get(), p.get(), M_inv.get());
768 gsl_set_error_handler(old_handler);
769
770 if (invert_status != GSL_SUCCESS) {
771 throw std::runtime_error("Matrix is singular");
772 }
773
774 return from_gsl_copy(M_inv.get());
775}
776
777RealMatrix RealMatrix::operator-() const {
778 auto new_data = std::vector(this->rows_ * this->cols_, 0.0);
779
780 for (size_t i = 0; i < this->rows_; ++i) {
781 for (size_t j = 0; j < this->cols_; ++j) {
782 new_data[i * cols_ + j] = -this->unchecked_at(i, j);
783 }
784 }
785
786 return RealMatrix(std::move(new_data), rows_, cols_);
787}
788
789RealMatrix &RealMatrix::operator+=(const RealMatrix &rhs) {
790 if (this->rows_ != rhs.rows() || this->cols_ != rhs.cols())
791 throw std::invalid_argument("Matrices should have the same shape.");
792
793 for (size_t i = 0; i < this->rows_; ++i) {
794 for (size_t j = 0; j < this->cols_; ++j) {
795 this->unchecked_at(i, j) += rhs.unchecked_at(i, j);
796 }
797 }
798
799 return *this;
800}
801
802RealMatrix &RealMatrix::operator-=(const RealMatrix &rhs) {
803 if (this->rows_ != rhs.rows() || this->cols_ != rhs.cols())
804 throw std::invalid_argument("Matrices should have the same shape.");
805
806 for (size_t i = 0; i < this->rows_; ++i) {
807 for (size_t j = 0; j < this->cols_; ++j) {
808 this->unchecked_at(i, j) -= rhs.unchecked_at(i, j);
809 }
810 }
811
812 return *this;
813}
814
815RealMatrix &RealMatrix::operator*=(const RealMatrix &rhs) {
816 if (this->cols_ != rhs.rows())
817 throw std::invalid_argument("Matrices don't have the right shape to be multiplied.");
818
819 const size_t M = rows_;
820 const size_t N = rhs.cols_;
821 const size_t K = cols_;
822
823 std::vector<double> result(M * N, 0.0);
824
825 constexpr size_t BLOCK = 32;
826
827 for (size_t ii = 0; ii < M; ii += BLOCK) {
828 for (size_t kk = 0; kk < K; kk += BLOCK) {
829 for (size_t jj = 0; jj < N; jj += BLOCK) {
830
831 size_t i_max = std::min(ii + BLOCK, M);
832 size_t k_max = std::min(kk + BLOCK, K);
833 size_t j_max = std::min(jj + BLOCK, N);
834
835 for (size_t i = ii; i < i_max; ++i) {
836 for (size_t k = kk; k < k_max; ++k) {
837 double aik = unchecked_at(i, k);
838 const double* rhs_row = &rhs.data[k * rhs.cols_];
839
840 double* res_row = &result[i * N];
841
842 for (size_t j = jj; j < j_max; ++j) {
843 res_row[j] += aik * rhs_row[j];
844 }
845 }
846 }
847 }
848 }
849 }
850
851 data = std::move(result);
852 cols_ = N;
853 return *this;
854}
855
856RealMatrix &RealMatrix::operator*=(double scalar) {
857 for (double& x: data)
858 x *= scalar;
859
860 return *this;
861}
862
863RealMatrix &RealMatrix::operator/=(double scalar) {
864 if (fpeq(scalar, 0.0))
865 throw std::invalid_argument("Division by zero.");
866
867 for (double& x: data)
868 x /= scalar;
869
870 return *this;
871}
872
873double &RealMatrix::unchecked_at(size_t i, size_t j) {
874 return data[i * cols_ + j];
875}
876
877const double &RealMatrix::unchecked_at(size_t i, size_t j) const {
878 return data[i * cols_ + j];
879}
nearest_psd
RealMatrix nearest_psd(RealMatrix R, double thr)
Projects a matrix to the nearest positive semi-definite correlation-like matrix.
Definition Matrix.cpp:475
eye
RealMatrix eye(std::size_t n)
Returns the identity matrix of size n.
Definition Matrix.cpp:453
diag
RealMatrix diag(const gsl_vector *X)
Builds a diagonal matrix from a GSL vector.
Definition Matrix.cpp:464
cholesky_L
RealMatrix cholesky_L(RealMatrix R)
Returns the lower-triangular Cholesky factor of a PSD matrix.
Definition Matrix.cpp:508
block_diag
RealMatrix block_diag(const std::vector< RealMatrix > &blocks)
Builds a block-diagonal matrix from several square blocks.
Definition Matrix.cpp:529
operator<<
std::ostream & operator<<(std::ostream &os, RealMatrix A)
Definition Matrix.cpp:439
Matrix.h
Lightweight dense real-matrix utilities built on top of STL storage and GSL backends.
make_gsl_vector
gsl_vector_sptr make_gsl_vector(size_t size)
Allocates a GSL vector with automatic RAII destruction.
Definition Matrix.h:79
make_gsl_permutation
gsl_permutation_sptr make_gsl_permutation(size_t size)
Allocates a GSL permutation with automatic RAII destruction.
Definition Matrix.h:88
gsl_matrix_sptr
std::unique_ptr< gsl_matrix, decltype(&gsl_matrix_free)> gsl_matrix_sptr
Owning smart pointer for a GSL matrix.
Definition Matrix.h:53
gsl_vector_sptr
std::unique_ptr< gsl_vector, decltype(&gsl_vector_free)> gsl_vector_sptr
Owning smart pointer for a GSL vector.
Definition Matrix.h:56
make_eigen_workspace
gsl_eigen_workspace_sptr make_eigen_workspace(size_t n)
Allocates a GSL symmetric eigensolver workspace with RAII destruction.
Definition Matrix.h:97
make_gsl_matrix
gsl_matrix_sptr make_gsl_matrix(size_t rows, size_t cols)
Allocates a GSL matrix with automatic RAII destruction.
Definition Matrix.h:70
gsl_eigen_workspace_sptr
std::unique_ptr< gsl_eigen_symmv_workspace, decltype(&gsl_eigen_symmv_free)> gsl_eigen_workspace_sptr
Owning smart pointer for a GSL symmetric-eigensystem workspace.
Definition Matrix.h:62
gsl_permutation_sptr
std::unique_ptr< gsl_permutation, decltype(&gsl_permutation_free)> gsl_permutation_sptr
Owning smart pointer for a GSL permutation.
Definition Matrix.h:59
RealMatrix::remove_row_and_column
void remove_row_and_column(std::size_t dim_idx)
Removes one row and one column with the same index.
Definition Matrix.cpp:639
RealMatrix::remove_column
void remove_column(std::size_t col_idx)
Removes one column from the matrix.
Definition Matrix.cpp:624
RealMatrix::operator/=
RealMatrix & operator/=(double scalar)
In-place scalar division.
Definition Matrix.cpp:863
RealMatrix::rows
std::size_t rows() const
Returns the number of rows.
Definition Matrix.cpp:601
RealMatrix::eig
EigenSystem eig() const
Computes the eigensystem of a symmetric matrix.
Definition Matrix.cpp:710
RealMatrix::at
double & at(size_t i, size_t j)
Returns a mutable reference to element (i,j) with bounds checking.
Definition Matrix.cpp:587
RealMatrix::operator-=
RealMatrix & operator-=(const RealMatrix &rhs)
In-place matrix subtraction.
Definition Matrix.cpp:802
RealMatrix::operator*=
RealMatrix & operator*=(const RealMatrix &rhs)
In-place matrix multiplication.
Definition Matrix.cpp:815
RealMatrix::RealMatrix
RealMatrix()
Constructs an empty 0×0 matrix.
Definition Matrix.h:149
RealMatrix::from_gsl_copy
static RealMatrix from_gsl_copy(const gsl_matrix *A)
Creates a RealMatrix by copying data from a GSL matrix.
Definition Matrix.cpp:658
RealMatrix::remove_row
void remove_row(std::size_t row_idx)
Removes one row from the matrix.
Definition Matrix.cpp:609
RealMatrix::slogdet
SignedLogDet slogdet() const
Computes the signed logarithmic determinant via LU decomposition.
Definition Matrix.cpp:731
RealMatrix::to_gsl_matrix
gsl_matrix_sptr to_gsl_matrix() const
Converts this matrix to a newly allocated GSL matrix.
Definition Matrix.cpp:672
RealMatrix::operator+=
RealMatrix & operator+=(const RealMatrix &rhs)
In-place matrix addition.
Definition Matrix.cpp:789
RealMatrix::unchecked_at
double & unchecked_at(size_t i, size_t j)
Returns a mutable reference to element (i,j) without bounds checking.
Definition Matrix.cpp:873
RealMatrix::is_symmetric
bool is_symmetric() const
Checks whether the matrix is symmetric.
Definition Matrix.cpp:684
RealMatrix::operator-
RealMatrix operator-() const
Unary minus.
Definition Matrix.cpp:777
RealMatrix::cols
std::size_t cols() const
Returns the number of columns.
Definition Matrix.cpp:605
RealMatrix::transpose
RealMatrix transpose() const
Returns the transpose of the matrix.
Definition Matrix.cpp:698
RealMatrix::inv
RealMatrix inv() const
Computes the inverse of the matrix via LU decomposition.
Definition Matrix.cpp:750
fpeq
std::enable_if_t< not std::numeric_limits< T >::is_integer, bool > fpeq(T, T, std::size_t n=10)
Compares two floating point numbers with a given precision.
Definition special_generic.h:107
EigenSystem
Container for an eigendecomposition.
Definition Matrix.h:489
EigenSystem::D
RealMatrix D
Definition Matrix.h:490
EigenSystem::P
RealMatrix P
Diagonal matrix of eigenvalues.
Definition Matrix.h:491
SignedLogDet
Signed logarithmic determinant representation.
Definition Matrix.h:117
SignedLogDet::sign
int sign
Natural logarithm of the absolute determinant.
Definition Matrix.h:119
SignedLogDet::logdet
double logdet
Definition Matrix.h:118