libcrn  3.9.5
A document image processing library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRNSquareMatrixDouble.cpp
Go to the documentation of this file.
1 /* Copyright 2008-2016 INSA-Lyon, ENS-Lyon
2  *
3  * This file is part of libcrn.
4  *
5  * libcrn is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU Lesser General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * libcrn is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public License
16  * along with libcrn. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * file: CRNSquareMatrixDouble.cpp
19  *
20  * \author Jean DUONG, Yann LEYDIER
21  */
22 
23 #include <CRNException.h>
25 #include <CRNMath/CRNMath.h>
27 #include <CRNData/CRNDataFactory.h>
28 
29 #include <stdio.h>
30 #include <iostream>
31 #include <limits>
32 #include <CRNi18n.h>
33 
34 using namespace crn;
35 
37  MatrixDouble(m)
38 {
39  if (rows != cols)
40  throw ExceptionDimension(StringUTF8("SquareMatrixDouble::SquareMatrixDouble(const Matrix<double> &m): ") + _("the matrix is not square."));
41 }
42 
44  MatrixDouble(std::move(m))
45 {
46  if (rows != cols)
47  throw ExceptionDimension(StringUTF8("SquareMatrixDouble::SquareMatrixDouble(const Matrix<double> &m): ") + _("the matrix is not square."));
48 }
49 
50 SquareMatrixDouble::SquareMatrixDouble(const std::vector<std::vector<double>> &m):
51  MatrixDouble(m)
52 {
53  if (rows != cols)
54  throw ExceptionDimension(StringUTF8("SquareMatrixDouble::SquareMatrixDouble(const std::vector<std::vector<double>> &m): ") + _("the matrix is not square."));
55 }
56 
57 SquareMatrixDouble::SquareMatrixDouble(std::vector<std::vector<double>> &&m):
58  MatrixDouble(std::move(m))
59 {
60  if (rows != cols)
61  throw ExceptionDimension(StringUTF8("SquareMatrixDouble::SquareMatrixDouble(const std::vector<std::vector<double>> &m): ") + _("the matrix is not square."));
62 }
63 
64 // Special matrix constructors
65 #define MAX_GAUSS_W 40
66 #define MULT 1
67 
76 {
77  if (sigma == 0)
78  return SquareMatrixDouble(1, 1);
79 
80  size_t d = MAX_GAUSS_W;
81  // compute the width of the matrix
82  for (size_t tmp = 1; tmp < MAX_GAUSS_W; ++tmp)
83  if (Gauss(double(tmp), sigma) < 0.1)
84  {
85  d = tmp;
86  break;
87  }
88  SquareMatrixDouble mat(d + d + 1);
89 
90  // fill central column
91  for (size_t r = 0; r < d; r++)
92  mat[r][d] = MULT * MeanGauss(double(d) - double(r), sigma);
93  for (size_t r = d; r < mat.rows; r++)
94  mat[r][d] = MULT * MeanGauss(double(r) - double(d), sigma);
95  // fill central row
96  for (size_t c = 0; c < d; c++)
97  mat[d][c] = MULT * MeanGauss(double(d) - double(c), sigma);
98  for (size_t c = d; c < mat.cols; c++)
99  mat[d][c] = MULT * MeanGauss(double(c) - double(d), sigma);
100  // fill first quadrant
101  for (size_t r = 0; r < d; r++)
102  for (size_t c = d + 1; c < mat.cols; c++)
103  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
104  // fill second quadrant
105  for (size_t r = 0; r < d; r++)
106  for (size_t c = 0; c < d; c++)
107  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
108  // fill third quadrant
109  for (size_t r = d + 1; r < mat.rows; r++)
110  for (size_t c = 0; c < d; c++)
111  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
112  // fill fourth quadrant
113  for (size_t r = d + 1; r < mat.rows; r++)
114  for (size_t c = d + 1; c < mat.cols; c++)
115  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
116 
117  return mat;
118 }
119 
128 {
129  if (sigma == 0)
130  {
131  SquareMatrixDouble m(3, 0);
132  m[1][0] = 1;
133  m[1][2] = -1;
134  return m;
135  }
136 
137  size_t d = MAX_GAUSS_W;
138  // compute the width of the matrix
139  for (size_t tmp = 1; tmp < MAX_GAUSS_W; tmp++)
140  if (Gauss(double(tmp), sigma) < 0.1)
141  {
142  d = tmp;
143  break;
144  }
145 
146  SquareMatrixDouble mat(d + d + 1);
147 
148  // fill central column
149  for (size_t r = 0; r < d; r++)
150  mat[r][d] = MULT * MeanGauss(double(d) - double(r), sigma);
151  for (size_t r = d; r < mat.rows; r++)
152  mat[r][d] = MULT * MeanGauss(double(r) - double(d), sigma);
153  // fill central row
154  for (size_t c = 0; c < d; c++)
155  mat[d][c] = MULT * MeanGauss(double(d) - double(c), sigma);
156  for (size_t c = d; c < mat.cols; c++)
157  mat[d][c] = MULT * MeanGauss(double(c) - double(d), sigma);
158  // fill first quadrant
159  for (size_t r = 0; r < d; r++)
160  for (size_t c = d + 1; c < mat.cols; c++)
161  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
162  // fill second quadrant
163  for (size_t r = 0; r < d; r++)
164  for (size_t c = 0; c < d; c++)
165  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
166  // fill third quadrant
167  for (size_t r = d + 1; r < mat.rows; r++)
168  for (size_t c = 0; c < d; c++)
169  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
170  // fill fourth quadrant
171  for (size_t r = d + 1; r < mat.rows; r++)
172  for (size_t c = d + 1; c < mat.cols; c++)
173  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
174 
175  // derivate
176  for (size_t r = 0; r < mat.rows; r++)
177  for (size_t c = 0; c < mat.cols; c++)
178  {
179  mat[r][c] *= double(d) - double(c);
180  mat[r][c] /= double(d);
181  }
182 
183  return mat;
184 }
185 
194 {
195  if (sigma == 0)
196  {
197  SquareMatrixDouble m(3, 0);
198  m[0][1] = 1;
199  m[2][1] = -1;
200  return m;
201  }
202 
203  size_t d = MAX_GAUSS_W;
204  // compute the width of the matrix
205  for (size_t tmp = 1; tmp < MAX_GAUSS_W; tmp++)
206  if (Gauss(double(tmp), sigma) < 0.1)
207  {
208  d = tmp;
209  break;
210  }
211  SquareMatrixDouble mat(d + d + 1);
212 
213  // fill central column
214  for (size_t r = 0; r < d; r++)
215  mat[r][d] = MULT * MeanGauss(double(d) - double(r), sigma);
216  for (size_t r = d; r < mat.rows; r++)
217  mat[r][d] = MULT * MeanGauss(double(r) - double(d), sigma);
218  // fill central row
219  for (size_t c = 0; c < d; c++)
220  mat[d][c] = MULT * MeanGauss(double(d) - double(c), sigma);
221  for (size_t c = d; c < mat.cols; c++)
222  mat[d][c] = MULT * MeanGauss(double(c) - double(d), sigma);
223  // fill first quadrant
224  for (size_t r = 0; r < d; r++)
225  for (size_t c = d + 1; c < mat.cols; c++)
226  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
227  // fill second quadrant
228  for (size_t r = 0; r < d; r++)
229  for (size_t c = 0; c < d; c++)
230  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
231  // fill third quadrant
232  for (size_t r = d + 1; r < mat.rows; r++)
233  for (size_t c = 0; c < d; c++)
234  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
235  // fill fourth quadrant
236  for (size_t r = d + 1; r < mat.rows; r++)
237  for (size_t c = d + 1; c < mat.cols; c++)
238  mat[r][c] = mat[r][d] * mat[d][c] / MULT;
239 
240  // derivate
241  for (size_t r = 0; r < mat.rows; r++)
242  for (size_t c = 0; c < mat.cols; c++)
243  {
244  mat[r][c] *= double(d) - double(r);
245  mat[r][c] /= double(d);
246  }
247 
248  return mat;
249 }
250 
257 {
258  bool Status = true;
259  size_t r = 1;
260 
261  while ((Status) && (r < rows))
262  {
263  size_t c = 0;
264 
265  while ((Status) && (c < r))
266  {
267  Status = (At(r, c) == 0.0);
268  ++c;
269  }
270 
271  ++r;
272  }
273 
274  return Status;
275 }
276 
283 {
284  bool Status = true;
285  size_t r = 0;
286 
287  while ((Status) && (r < rows - 1))
288  {
289  size_t c = r + 1;
290 
291  while ((Status) && (c < cols))
292  {
293  Status = (At(r, c) == 0.0);
294  ++c;
295  }
296 
297  ++r;
298  }
299 
300  return Status;
301 }
302 
309 {
310  auto status = true;
311  auto c = 0;
312 
313  while (status && (c < int(cols) - 2))
314  {
315  auto r = c + 2;
316 
317  while (status && (r < int(rows)))
318  {
319  status = status && !At(r, c);
320  ++r;
321  }
322 
323  ++c;
324  }
325 
326  return status;
327 }
328 
335 {
336  bool Status = true;
337  size_t r = 0;
338 
339  while ((Status) && (r < rows - 1))
340  {
341  size_t c = r + 1;
342 
343  while ((Status) && (c < cols))
344  {
345  Status = (At(r, c) == 0.0) && (At(c, r) == 0.0);
346  ++c;
347  }
348 
349  ++r;
350  }
351 
352  return Status;
353 }
354 
361 {
362  double sum = 0;
363 
364  for (size_t r = 0; r < rows; ++r)
365  sum += At(r, r);
366 
367  return sum;
368 }
369 
376 {
377  double P = 1.0;
378 
379  for (size_t r = 0; r < rows; ++r)
380  P *= At(r, r);
381 
382  return P;
383 }
384 
389 {
390  for (size_t l = 0; l < rows - 1; ++l)
391  for (size_t c = l + 1; c < cols; ++c)
392  std::swap(At(l, c), At(c, l));
393  return *this;
394 }
395 
407 {
408  if ((r >= rows) || (c >= cols))
409  throw ExceptionDomain(StringUTF8("SquareMatrixDouble::MakeMinor(size_t r, size_t c): ") +
410  _("row or column index out of range"));
411 
412  SquareMatrixDouble M(rows - 1);
413 
414  size_t i = 0;
415 
416  for (size_t rindex = 0; rindex < rows; ++rindex)
417  {
418  if (rindex != r)
419  {
420  size_t j = 0;
421 
422  for (size_t cindex = 0; cindex < cols; ++cindex)
423  {
424  if (cindex != c)
425  {
426  M[i][j] = At(rindex, cindex);
427  ++j;
428  }
429  }
430 
431  ++i;
432  }
433  }
434 
435  return M;
436 }
437 
448 double SquareMatrixDouble::Cofactor(size_t r, size_t c) const
449 {
450  if ((r >= rows) || (c >= cols))
451  throw ExceptionDomain(StringUTF8("SquareMatrixDouble::Cofactor(size_t r, size_t c): ") +
452  _("row or column index out of range"));
453 
454  SquareMatrixDouble M = MakeMinor(r, c);
455  double d = M.Determinant();
456 
457  if ((r + c) % 2 == 1) // Parity test
458  d = -d;
459 
460  return d;
461 }
462 
469 {
470  // Case 1x1 matrix
471  if (rows == 1)
472  return At(0, 0);
473 
474  // Case 2x2 matrix
475  if (rows == 2)
476  return At(0, 0) * At(1, 1) - At(0, 1) * At(1, 0);
477 
478  // Case 3x3 matrix (Sarrus rule)
479  if (rows == 3)
480  {
481  double det = 0.0;
482 
483  det += At(0, 0) * At(1, 1) * At(2, 2);
484  det += At(1, 0) * At(2, 1) * At(0, 2);
485  det += At(2, 0) * At(0, 1) * At(1, 2);
486  det -= At(0, 2) * At(1, 1) * At(2, 0);
487  det -= At(1, 2) * At(2, 1) * At(0, 0);
488  det -= At(2, 2) * At(0, 1) * At(1, 0);
489 
490  return det;
491  }
492 
493  // Case triangular matrices
495  return DiagonalProduct();
496 
498  // General case //
500 
501  double det = 0.0;
502 
503  size_t r = 0;
504  size_t c = 0;
505  size_t MaxNullInRow = 0;
506  size_t MaxNullInColumn = 0;
507 
508  // Search most sparse row (r) and most sparse column (c)
509 
510  for (size_t k = 0; k < rows; ++k)
511  {
512  size_t nr = CountNullCellsInRow(k);
513  size_t nc = CountNullCellsInColumn(k);
514 
515  if (nr > MaxNullInRow)
516  {
517  MaxNullInRow = nr;
518  r = k;
519  }
520 
521  if (nc > MaxNullInColumn)
522  {
523  MaxNullInColumn = nc;
524  c = k;
525  }
526  }
527 
528  // Expansion by minors
529 
530  if (MaxNullInRow > MaxNullInColumn)
531  {
532  for (size_t k = 0; k < cols; ++k)
533  {
534  double Erk = At(r, k);
535 
536  if (Erk != 0.0)
537  det += Erk * Cofactor(r, k);
538  }
539  }
540  else
541  {
542  for (size_t k = 0; k < rows; k++)
543  {
544  double Ekc = At(k, c);
545 
546  if (Ekc != 0.0)
547  det += Ekc * Cofactor(k, c);
548  }
549  }
550 
551  return det;
552 }
553 
560 {
561  double dt = Determinant();
562 
563  if (dt != 0.0)
564  {
566 
567  for (size_t r = 0; r < rows; ++r)
568  for (size_t c = 0; c < cols; ++c)
569  C[r][c] = Cofactor(r, c);
570 
571  C.Transpose();
572  C *= 1.0 / dt;
573 
574  return C;
575  }
576 
577  throw ExceptionRuntime(_("The matrix cannot be inversed."));
578 }
579 
588 {
589  SquareMatrixDouble id(n);
590 
591  for (size_t r = 0; r < n; ++r)
592  for (size_t c = 0; c < n; ++c)
593  {
594  if (r == c)
595  id[r][c] = 1.0;
596  else
597  id[r][c] = 0.0;
598  }
599 
600  return id;
601 }
602 
611 {
612  size_t n = rows;
614  SquareMatrixDouble M = *this;
615 
616  // Transform this matrix to triangular matrix
617 
618  for (size_t c = 0; c < n - 1; ++c)
619  {
620  // Search the greatest pivot in column
621 
622  double Pivot = M[c][c];
623  double AbsMaxPivot = fabs(Pivot);
624  size_t RowIndex = c;
625 
626  for (size_t r = c + 1 ; r < n; ++r)
627  {
628  double Candidate = M[r][c];
629 
630  if (fabs(Candidate) > AbsMaxPivot)
631  {
632  Pivot = Candidate;
633  AbsMaxPivot = fabs(Pivot);
634  RowIndex = r;
635  }
636  }
637 
638  // If no non-null pivot found, matrix is non inversible
639 
640  if (Pivot == 0.0)
641  throw ExceptionRuntime(crn::StringUTF8(_("Matrix cannot be inversed. Null pivot at column: ")) + c);
642 
643  if (RowIndex != c)
644  {
645  M.SwapRows(c, RowIndex);
646  Id.SwapRows(c, RowIndex);
647  }
648 
649  // Elimination
650 
651  for (size_t r = c + 1; r < n; ++r)
652  {
653  double Coeff = M[r][c];
654 
655  if (Coeff != 0.0)
656  {
657  double Scale = - Coeff / Pivot;
658 
659  for (size_t k = 0; k < n; ++k)
660  {
661  M[r][k] += M[c][k] * Scale;
662  Id[r][k] += Id[c][k] * Scale;
663  }
664  }
665  }
666  }
667 
668  // Transform this matrix to diagonal
669 
670  for (int c = int(n) - 1; c > 0; --c)
671  {
672  double DCoeff = M[c][c];
673 
674  if (DCoeff != 0.0)
675  {
676  for (int r = c - 1; r >= 0; --r)
677  {
678  double Coeff = M[r][c];
679 
680  if (Coeff != 0.0)
681  {
682  double Scale = - Coeff / DCoeff;
683 
684  for (size_t k = 0; k < n; ++k)
685  {
686  Id[r][k] += Id[c][k] * Scale;
687  M[r][k] += M[c][k] * Scale;
688  }
689  }
690  }
691  }
692  }
693 
694  // Transform to id
695 
696  for (size_t r = 0; r < n; r++)
697  Id.MultRow(r, 1.0 / M[r][r]);
698 
699  return Id;
700 }
701 
714 {
715  size_t n = GetRows();
716  SquareMatrixDouble L(n, 0.0);
717 
718  L[0][0] = sqrt(At(0, 0));
719 
720  if (n > 2)
721  {
722  for (int j = 0; j < int(n); ++j)
723  {
724  auto Sum = 0.0;
725 
726  for (int k = 0; k <= j - 1; ++k)
727  Sum += Sqr(L[j][k]);
728 
729  L[j][j] = sqrt(At(j, j) - Sum);
730 
731  for (int i = j + 1 ; i < int(n); ++i)
732  {
733  Sum = 0.0;
734 
735  for (int k = 0; k <= j - 1; ++k)
736  Sum += L[i][k] * L[j][k];
737 
738  L[i][j] = (At(i, j) - Sum) / L[j][j];
739  }
740  }
741  }
742  else
743  {
744  L[1][0] = At(0, 1) / L[0][0];
745  L[1][1] = sqrt(At(1, 1) - Sqr(At(0, 1)) / At(0, 0));
746  }
747 
748  return L;
749 }
750 
757 {
758  const auto n = Min(int(rows), int(cols)); // Vector space dimension
759  auto hess_mat(*this);
760 
761  std::vector<std::vector<double>> reflector_seeds;
762 
763  reflector_seeds.reserve(n - 2);
764 
765  for (auto k = 0; k <= n - 3; ++k)
766  {
767  // Generate the Househoder reflector Pk
768  std::vector<double> u_k;
769 
770  u_k.reserve(n - k);
771 
772  for (auto i = k + 1; i < n; ++i)
773  u_k.push_back(hess_mat.At(i, k));
774 
775  u_k[0] += SignOf(u_k[0]) * Pythagoras(u_k.begin(), u_k.end());
776  Scale(u_k.begin(), u_k.end(), 1.0 / Pythagoras(u_k.begin(), u_k.end()));
777  reflector_seeds.push_back(u_k);
778 
779  // Left product by P_k
780  for (auto j = k; j < n; ++j)
781  {
782  auto scale_current_column = 0.0;
783  auto it = u_k.begin();
784 
785  for (auto i = k + 1; i < n; ++i)
786  {
787  scale_current_column += (*it) * hess_mat.At(i, j);
788  ++it;
789  }
790 
791  it = u_k.begin();
792 
793  for (auto i = k + 1; i < n; ++i)
794  {
795  hess_mat.At(i, j) -= 2 * (*it) * scale_current_column;
796  ++it;
797  }
798  }
799 
800  // Right product by P_k
801  for (auto i = 0; i < n; ++i)
802  {
803  auto scale_current_row = 0.0;
804  auto it = u_k.begin();
805 
806  for (auto j = k + 1; j < n; ++j)
807  {
808  scale_current_row += hess_mat.At(i, j) * (*it);
809  ++it;
810  }
811 
812  it = u_k.begin();
813 
814  for (auto j = k + 1; j < n; ++j)
815  {
816  hess_mat.At(i, j) -= 2 * scale_current_row * (*it);
817  ++it;
818  }
819  }
820  }
821 
822  return hess_mat;
823 }
824 
825 
834 std::multimap<double, MatrixDouble> SquareMatrixDouble::MakeSpectralEigensystem() const
835 {
836  std::multimap<double, MatrixDouble> eigen_pairs;
837 
838  auto g_1 = At(0, 0);
839  auto g_2 = At(0, 1);
840  auto g_3 = At(1, 1);
841 
842  MatrixDouble eigen_vector_1((size_t)2, (size_t)1);
843  MatrixDouble eigen_vector_2((size_t)2, (size_t)1);
844 
845  if (g_2 == 0.0) // Trivial case (matrix already diagonal)
846  {
847  eigen_vector_1[0][0] = 1.0;
848  eigen_vector_1[1][0] = 0;
849  eigen_vector_2[0][0] = 0.0;
850  eigen_vector_2[1][0] = 1.0;
851 
852  eigen_pairs.insert(std::make_pair(g_1, eigen_vector_1));
853  eigen_pairs.insert(std::make_pair(g_3, eigen_vector_2));
854  }
855  else
856  {
857  auto delta = Sqr(g_1 - g_3) + 4 * Sqr(g_2);
858  auto sqrt_delta = sqrt(delta);
859 
860  auto eigen_value_1 = (g_1 + g_3 + sqrt_delta) / 2.0;
861  auto eigen_value_2 = (g_1 + g_3 - sqrt_delta) / 2.0;
862 
863  eigen_vector_1[0][0] = 1.0;
864  eigen_vector_1[1][0] = (eigen_value_1 - g_1) / g_2;
865 
866  // For symetric matrix, eigenvectors are orthogonal: (u, v) and (-v, u)
867 
868  eigen_vector_2[0][0] = - eigen_vector_1[1][0];
869  eigen_vector_2[1][0] = eigen_vector_1[0][0];
870 
871  eigen_pairs.insert(std::make_pair(eigen_value_1, eigen_vector_1));
872  eigen_pairs.insert(std::make_pair(eigen_value_2, eigen_vector_2));
873  }
874 
875  return eigen_pairs;
876 }
877 
893 std::multimap<double, MatrixDouble> SquareMatrixDouble::MakeJacobiEigensystem(size_t MaxIteration) const
894 {
895  size_t MaxIter = 100;
896 
897  size_t i, j, k, n = GetRows();
898 
899  SquareMatrixDouble Eigenvectors(n, 0.0);
900  SquareMatrixDouble A = *this;
901 
902  for (k = 0; k < n; k++)
903  Eigenvectors[k][k] = 1.0;
904 
905  k = 0;
906  bool Loop = true;
907 
908  // Some variable matrices.
909 
910  SquareMatrixDouble Omega(n, 0.0);
911  SquareMatrixDouble NewA(n, 0.0);
912  SquareMatrixDouble B(n, 0.0);
913 
915  // Loop //
917 
918  while ((k < MaxIter) && Loop)
919  {
920  size_t p = 1, q = 0;
921 
922  // Search optimal p and q indices.
923 
924  double M = 0.0;
925 
926  for (i = 0; i < n; ++i)
927  for (j = 0; j < n; ++j)
928  if (i != j)
929  {
930  double a = fabs(A[i][j]);
931 
932  if (a > M)
933  {
934  M = a;
935  p = i;
936  q = j;
937  }
938  }
939 
940  // Build rotation.
941 
942  Omega.SetAll(0.0);
943 
944  for (i = 0; i < n; i++)
945  Omega[i][i] = 1.0;
946 
947  double Xi = (A[q][q] - A[p][p]) / (2.0 * A[p][q]);
948  double t = 1.0;
949 
950  if (Xi != 0.0)
951  {
952  std::set<double> Roots = QuadraticEquation::RealRoots(1.0, 2.0 * Xi, -1.0);
953 
954  t = Min(*Roots.begin(), *Roots.rbegin());
955  }
956 
957  double c = 1.0 / sqrt(1 + t * t);
958  double s = t / sqrt(1 + t * t);
959 
960  Omega[p][p] = c;
961  Omega[q][q] = c;
962  Omega[p][q] = s;
963  Omega[q][p] = -s;
964 
965  // Iterate eigenvectors.
966 
967  Eigenvectors *= Omega;
968 
969  // Iterate A.
970 
971  B = A * Omega;
972  Omega.Transpose();
973  NewA = Omega * B;
974  Loop = (NewA != A);
975  A = NewA;
976  ++k;
977  }
978 
980  // End of main loop. //
982 
983  // Pack pairs of eigenvalues with corresponding eigenvectors
984  // in descending order according to eigenvalues.
985 
986  std::multimap<double, MatrixDouble> EigenPairs;
987  std::vector<bool> Flags(n);
988 
989  for (k = 0; k < n; ++k)
990  Flags[k] = true;
991 
992  for (i = 0; i < n; ++i)
993  {
994  double MaxEigenvalue = -1.0; // Eigenvalues are all positive
995  size_t Argmax = 0;
996 
997  for (j = 0; j < n; ++j)
998  if (Flags[j])
999  {
1000  double Value = A[j][j];
1001 
1002  if (Value > MaxEigenvalue)
1003  {
1004  MaxEigenvalue = Value;
1005  Argmax = j;
1006  }
1007  }
1008 
1009  Flags[Argmax] = false;
1010 
1011  EigenPairs.insert(std::make_pair(A[Argmax][Argmax], Eigenvectors.MakeColumn(Argmax)));
1012  }
1013 
1014  return EigenPairs;
1015 }
1016 
1023 std::multimap<double, MatrixDouble> SquareMatrixDouble::MakeTQLIEigensystem(size_t maxiter) const
1024 {
1025  SquareMatrixDouble z(1);
1026  std::vector<double> eigenvalues, offdiag;
1027  tred2(z, eigenvalues, offdiag);
1028 
1029  // tqli
1030  for (int i = 1; i < int(rows); ++i) // renumber the off-diagonal
1031  offdiag[i - 1] = offdiag[i];
1032  offdiag[rows - 1] = 0.0;
1033  for (int l = 0; l < int(rows); ++l)
1034  {
1035  size_t iter = 0;
1036  int m;
1037  do
1038  {
1039  for (m = l; m < int(rows) - 1; ++m)
1040  { // look for a single small subdiagonal element to split the matrix
1041  double dd = Abs(eigenvalues[m]) + Abs(eigenvalues[m + 1]);
1042  if (Abs(offdiag[m]) <= std::numeric_limits<double>::epsilon() * dd)
1043  break;
1044  }
1045  if (m != l)
1046  {
1047  if (iter++ >= maxiter)
1048  throw ExceptionRuntime(StringUTF8("std::multimap<double, SMatrixDouble> SquareMatrixDouble::MakeTQLIEigensystem(size_t maxiter) const: ") + _("Too many iterations."));
1049  double g = (eigenvalues[l + 1] - eigenvalues[l]) / (2.0 * offdiag[l]);
1050  double r = Pythagoras(g, 1.0);
1051  g = eigenvalues[m] - eigenvalues[l] + offdiag[l] / (g + ((g < 0) ? -Abs(r) : Abs(r)));
1052  double s = 1.0, c = 1.0, p = 0.0;
1053  int i;
1054  for (i = m - 1; i >= l; --i)
1055  { // plane rotation followed by Givens rotations to restore tridiagonal form
1056  double f = s * offdiag[i];
1057  double b = c * offdiag[i];
1058  r = Pythagoras(f, g);
1059  offdiag[i + 1] = r;
1060  if (r == 0.0)
1061  { // recover from underflow
1062  eigenvalues[i + 1] -= p;
1063  offdiag[m] = 0.0;
1064  break;
1065  }
1066  s = f / r;
1067  c = g / r;
1068  g = eigenvalues[i + 1] - p;
1069  r = (eigenvalues[i] - g) * s + 2.0 * c * b;
1070  p = s * r;
1071  eigenvalues[i + 1] = g + p;
1072  g = c * r - b;
1073  for (int k = 0; k < int(rows); ++k)
1074  {
1075  f = z[k][i + 1];
1076  z[k][i + 1] = s * z.At(k, i) + c * f;
1077  z.At(k, i) = c * z.At(k, i) - s * f;
1078  }
1079  }
1080  if ((r == 0.0) && (i >= l))
1081  continue;
1082  eigenvalues[l] -= p;
1083  offdiag[l] = g;
1084  offdiag[m] = 0.0;
1085  }
1086  } while (m != l);
1087  }
1088 
1089  // export results
1090  std::multimap<double, MatrixDouble> out;
1091  for (size_t tmp = 0; tmp < eigenvalues.size(); ++tmp)
1092  {
1093  MatrixDouble ev(rows, 1, 0.0);
1094  for (size_t r = 0; r < rows; ++r)
1095  ev[r][0] = z[r][tmp];
1096  out.insert(std::make_pair(eigenvalues[tmp], ev));
1097  }
1098  return out;
1099 }
1100 
1106 void SquareMatrixDouble::tred2(SquareMatrixDouble &z, std::vector<double> &diag, std::vector<double> &offdiag) const
1107 {
1108  // initialize output variables
1109  z = *this;
1110  std::vector<double>(rows).swap(diag);
1111  std::vector<double>(rows).swap(offdiag);
1112 
1113  for (int i = int(rows) - 1; i > 0; --i)
1114  {
1115  int l = i - 1;
1116  double h = 0.0;
1117  if (l > 0)
1118  {
1119  double scale = 0.0;
1120  for (int k = 0; k < i; ++k)
1121  scale += Abs(z.At(i, k));
1122  if (scale == 0.0) // skip transformation
1123  offdiag[i] = z.At(i, l);
1124  else
1125  {
1126  for (int k = 0; k < i; ++k)
1127  {
1128  z.At(i, k) /= scale;
1129  h += Sqr(z.At(i, k));
1130  }
1131  double f = z.At(i, l);
1132  double g = (f > 0.0) ? -sqrt(h) : sqrt(h);
1133  offdiag[i] = scale * g;
1134  h -= f * g;
1135  z.At(i, l) = f - g;
1136  f = 0.0;
1137  for (int j = 0; j < i; ++j)
1138  {
1139  z.At(j, i) = z.At(i, j) / h;
1140  g = 0.0;
1141  for (int k = 0; k < j + 1; ++k)
1142  g += z.At(j, k) * z.At(i, k);
1143  for (int k = j + 1; k < i; ++k)
1144  g += z.At(k, j) * z.At(i, k);
1145  offdiag[j] = g / h;
1146  f += offdiag[j] * z.At(i, j);
1147  }
1148  double hh = f / (h + h);
1149  for (int j = 0; j < i; ++j)
1150  {
1151  f = z.At(i, j);
1152  offdiag[j] = g = offdiag[j] - hh * f;
1153  for (int k = 0; k < j + 1; ++k)
1154  z.At(j, k) -= f * offdiag[k] + g * z.At(i, k);
1155  }
1156  }
1157  }
1158  else
1159  offdiag[i] = z.At(i, l);
1160  diag[i] = h;
1161  }
1162  diag[0] = 0.0;
1163  offdiag[0] = 0.0;
1164  for (size_t i = 0; i < rows; ++i)
1165  { // begin accumulation of transformation matrices
1166  if (diag[i] != 0.0)
1167  {
1168  for (size_t j = 0; j < i; ++j)
1169  {
1170  double g = 0.0;
1171  for (size_t k = 0; k < i; ++k)
1172  g += z.At(i, k) * z.At(k, j);
1173  for (size_t k = 0; k < i; ++k)
1174  z.At(k, j) -= g * z.At(k, i);
1175  }
1176  }
1177  diag[i] = z.At(i, i);
1178  z.At(i, i) = 1.0; // reset row and column of z to identity matrix for next iteration
1179  for (size_t j = 0; j < i; ++j)
1180  z.At(j, i) = z.At(i, j) = 0.0;
1181  }
1182 
1183 }
1184 
1185 
1195 std::vector<std::complex<double>> SquareMatrixDouble::Eigenvalues(size_t max_iter) const
1196 {
1197  const auto n = int(GetCols());
1198 
1199  if (IsUpperTriangular())
1200  {
1201  std::vector<std::complex<double>> eigenvalues;
1202 
1203  for (auto k = 0; k < n; ++k)
1204  eigenvalues.push_back(std::complex<double>(At(k, k)));
1205 
1206  return eigenvalues;
1207  }
1208 
1209  // Implementation done for [1..n][1..n] matrices.
1210  // Some tricky index shift needed for [0..n-1][0..n-1] matrices.
1211  // Still to be done!
1212  // For now we just store [0..n-1][0..n-1] input matrix in a wider
1213  // [0..n][0..n] matrix to comply with code proposed in
1214  // "Numerical Recipes in C" (0th row and 0th colums not touched)
1215  SquareMatrixDouble a(n + 1, 0.0);
1216 
1217  if (IsUpperHessenberg())
1218  for (auto r = 0; r < n; ++r)
1219  for (auto c = 0; c < n; ++c)
1220  a[r + 1][c + 1] = At(r, c);
1221  else
1222  {
1223  const auto& hess_mat = MakeUpperHessenberg();
1224 
1225  for (auto r = 0; r < n; ++r)
1226  for (auto c = 0; c < n; ++c)
1227  a[r + 1][c + 1] = hess_mat.At(r, c);
1228  }
1229 
1230  std::vector<double> wr(n + 1);
1231  std::vector<double> wi(n + 1);
1232 
1233  int nn = n, m, l;
1234  double z, y, x, w, v, u, t = 0.0, s, r, q, p;
1235 
1236  // Compute matrix norm for possible use in
1237  // locating single small subdiagonal element.
1238  auto anorm = 0.0;
1239 
1240  for (auto i = 1; i <= n; ++i)
1241  for (auto j = Max(i - 1, 1); j <= n; ++j)
1242  anorm += Abs(a[i][j]);
1243 
1244  // Gets changed only by an exceptional shift.
1245  while (nn >= 1)
1246  {
1247  // Begin search for next eigenvalue.
1248  int its = 0;
1249 
1250  do
1251  {
1252  for (l = nn; l >= 2; --l)
1253  { // Begin iteration: look for single small subdiagonal element.
1254  s = Abs(a[l - 1][l - 1]) + Abs(a[l][l]);
1255 
1256  if (s == 0.0)
1257  s = anorm;
1258  if (Abs(a[l][l - 1]) + s == s) // Because s and a[l][l-1] may have very different magnitudes
1259  break;
1260  }
1261 
1262  x = a[nn][nn];
1263 
1264  if (l == nn)
1265  { // One root found.
1266  wr[nn] = x + t;
1267  wi[nn--] = 0.0;
1268  }
1269  else
1270  {
1271  y = a[nn - 1][nn - 1];
1272  w = a[nn][nn - 1] * a[nn - 1][nn];
1273 
1274  if (l == (nn - 1))
1275  { // Two roots found...
1276  p = (y - x) / 2.0;
1277  q = Sqr(p) + w;
1278  z = sqrt(Abs(q));
1279  x += t;
1280 
1281  if (q >= 0.0)
1282  { // ... a real pair.
1283  //z = p + SIGN(z, p);
1284  z = p + Abs(z) * SignOf(p);
1285  wr[nn - 1] = wr[nn] = x + z;
1286 
1287  if (z)
1288  wr[nn] = x - w / z;
1289 
1290  wi[nn - 1] = wi[nn] = 0.0;
1291  }
1292  else
1293  { // ... a complex pair.
1294  wr[nn - 1] = wr[nn] = x + p;
1295  wi[nn - 1] = -(wi[nn] = z);
1296  }
1297 
1298  nn -= 2;
1299  }
1300  else
1301  { // No roots found. Continue iteration.
1302  if (its == 30)
1303  std::cout << "ARGH !!!" << std::endl;
1304 
1305  if (!(its % 10))
1306  { // Form exceptional shift.
1307  t += x;
1308 
1309  for (auto i = 1; i <= nn; ++i)
1310  a[i][i]-= x;
1311 
1312  s = Abs(a[nn][nn - 1]) + Abs(a[nn - 1][nn - 2]);
1313  y = x = 0.75 * s;
1314  w = -0.4375 * Sqr(s);
1315  }
1316 
1317  ++its;
1318 
1319  for (m = (nn - 2); m >= l; --m)
1320  { // Form shift and then look for 2 consecutive
1321  // small subdiagonal elements.
1322  z = a[m][m];
1323  r = x - z;
1324  s = y - z;
1325  p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
1326  q = a[m + 1][m + 1] - z - r - s;
1327  r = a[m + 2][m + 1];
1328  s = Abs(p) + Abs(q) + Abs(r);
1329 
1330  // Scale to prevent overflow or underflow.
1331  p /= s;
1332  q /= s;
1333  r /= s;
1334 
1335  if (m == l)
1336  break;
1337 
1338  u = Abs(a[m][m - 1]) * (Abs(q) + Abs(r));
1339  v = Abs(p) * (Abs(a[m - 1][m - 1]) + Abs(z) + Abs(a[m + 1][m + 1]));
1340 
1341  if (u + v == v)
1342  break;
1343  }
1344 
1345  for (auto i = m + 2; i <= nn; ++i)
1346  {
1347  a[i][i - 2] = 0.0;
1348 
1349  if (i != (m + 2))
1350  a[i][i - 3] = 0.0;
1351  }
1352 
1353  for (auto k = m; k <= nn - 1; ++k)
1354  { // Double QR step on rows l to nn
1355  // and columns m to nn
1356  if (k != m)
1357  {
1358  p = a[k][k - 1];
1359 
1360  // Begin setup of Householder vector
1361  q = a[k + 1][k - 1];
1362  r = 0.0;
1363 
1364  if (k != (nn - 1))
1365  r = a[k + 2][k - 1];
1366 
1367  // Scale to prevent overflow or underflow
1368  if ((x = Abs(p) + Abs(q) + Abs(r)) != 0.0)
1369  {
1370  p /= x;
1371  q /= x;
1372  r /= x;
1373  }
1374  }
1375 
1376  //if ((s = SIGN(Pythagoras(p, q, r) ,p)) != 0.0)
1377  if ((s = Pythagoras(p, q, r) * SignOf(p)) != 0.0)
1378  {
1379  if (k == m)
1380  {
1381  if (l != m)
1382  a[k][k - 1] = -a[k][k - 1];
1383  }
1384  else
1385  a[k][k - 1] = -s * x;
1386 
1387  p += s;
1388  x = p / s;
1389  y = q / s;
1390  z = r / s;
1391  q /= p;
1392  r /= p;
1393 
1394  for (auto j = k; j <= nn; ++j)
1395  { // Row modifcation.
1396  p = a[k][j] + q * a[k + 1][j];
1397 
1398  if (k != (nn - 1))
1399  {
1400  p += r * a[k + 2][j];
1401  a[k + 2][j] -= p * z;
1402  }
1403 
1404  a[k + 1][j] -= p * y;
1405  a[k][j] -= p * x;
1406  }
1407 
1408  auto mmin = nn < k + 3 ? nn : k + 3;
1409 
1410  for (auto i = l; i <= mmin; ++i)
1411  { // Column modifcation.
1412  p = x * a[i][k] + y * a[i][k + 1];
1413 
1414  if (k != (nn - 1))
1415  {
1416  p += z*a[i][k + 2];
1417  a[i][k + 2] -= p * r;
1418  }
1419 
1420  a[i][k + 1] -= p * q;
1421  a[i][k] -= p;
1422  }
1423  }
1424  }
1425  }
1426  }
1427  }
1428  while (l < nn - 1);
1429  }
1430 
1431  std::vector<std::complex<double>> eigenvalues;
1432 
1433  for (auto k = 1; k <=n; ++k)
1434  eigenvalues.push_back(std::complex<double>(wr[k], wi[k]));
1435 
1436  // Sort eigenvalues by magnitude
1437 
1438  std::sort(eigenvalues.begin(), eigenvalues.end(), [](const std::complex<double> &c1, const std::complex<double> &c2)
1439  {
1440  return norm(c1) < norm(c2);
1441  });
1442 
1443  return eigenvalues;
1444 }
1445 
1447  CRN_DATA_FACTORY_REGISTER(U"SquareMatrixDouble", SquareMatrixDouble)
1448  Cloner::Register<SquareMatrixDouble>();
1449 CRN_END_CLASS_CONSTRUCTOR(SquareMatrixDouble)
1450 
std::vector< std::complex< double > > Eigenvalues(size_t max_iter=30) const
Extract eigenvalues for matrix of real (eigenvalues may be complex)
static SquareMatrixDouble NewGaussianSobelX(double sigma)
Create Gaussian Sobel X derivation mask.
bool IsDiagonal() const
Check if matrix is diagonal.
bool IsUpperTriangular() const
Check if matrix is upper triangular.
size_t GetRows() const noexcept
Returns the number of rows.
Definition: CRNMatrix.h:157
size_t GetCols() const noexcept
Returns the number of columns.
Definition: CRNMatrix.h:163
SquareMatrixDouble MakeGaussJordanInverse() const
Invert using Gauss-Jordan elimination.
crn::StringUTF8 Id
Definition: CRNAltoUtils.h:31
virtual SquareMatrixDouble & Transpose() override
Transposition.
A generic runtime error.
Definition: CRNException.h:131
SquareMatrixDouble(size_t size, double value=0.0)
Constructor.
#define _(String)
Definition: CRNi18n.h:51
double MeanGauss(double x, double sigma)
Computes Gauss function at x for a given standard deviation (centered in 0) – to use with matrices...
Definition: CRNMath.h:81
const T & Max(const T &a, const T &b)
Returns the max of two values.
Definition: CRNMath.h:47
Matrix MakeColumn(size_t k) const
Extracts one column from matrix.
Definition: CRNMatrix.h:513
size_t CountNullCellsInColumn(size_t c) const
Counts null cells in column.
Definition: CRNMatrix.h:573
size_t CountNullCellsInRow(size_t r) const
Counts null cells in row.
Definition: CRNMatrix.h:561
SquareMatrixDouble MakeUpperHessenberg() const
Get the upper Hessenberg form of matrix.
#define CRN_END_CLASS_CONSTRUCTOR(classname)
Defines a class constructor.
Definition: CRNObject.h:198
double Trace() const
Trace.
#define MAX_GAUSS_W
std::multimap< double, MatrixDouble > MakeJacobiEigensystem(size_t MaxIteration=100) const
Perform diagonalization for symmetric matrix.
static SquareMatrixDouble NewGaussian(double sigma)
Create Gaussian matrix.
A generic domain error.
Definition: CRNException.h:83
void Scale(ITER it_begin, ITER it_end, const double s)
Scale a collection of numbers.
Definition: CRNMath.h:103
static SquareMatrixDouble NewIdentity(size_t n)
Create identity matrix.
void MultRow(size_t r, double v)
Scale one row from matrix.
Definition: CRNMatrix.h:320
bool IsUpperHessenberg() const
Check if matrix is upper Hessenberg.
SquareMatrixDouble MakeInverse() const
Invert using determinants.
std::multimap< double, MatrixDouble > MakeTQLIEigensystem(size_t maxiter=30) const
Perform diagonalization for positive symmetric matrix.
std::set< double > RealRoots(double a, double b, double c) noexcept
Real roots of trinom .
#define CRN_DATA_FACTORY_REGISTER(elemname, classname)
Registers a class to the data factory.
constexpr SumType< T > Sqr(const T &v) noexcept(noexcept(v *v))
Returns the square of a value.
Definition: CRNMath.h:61
void SetAll(const T &v)
Set all elements.
Definition: CRNMatrix.h:173
A dimension error.
Definition: CRNException.h:119
double matrix class
double Gauss(double x, double sigma)
Computes Gauss function at x for a given standard deviation (centered in 0)
Definition: CRNMath.h:78
constexpr int SignOf(const T &x) noexcept(noexcept(x< 0))
Returns the sign (-1 or 1) of a value.
Definition: CRNMath.h:59
void Abs(Image< T > &img, typename std::enable_if< std::is_arithmetic< T >::value >::type *dummy=nullptr) noexcept
Replaces each pixel by its absolute value.
Definition: CRNImageGray.h:47
const T & Min(const T &a, const T &b)
Returns the min of two values.
Definition: CRNMath.h:49
const double & At(size_t pos) const noexcept
Definition: CRNMatrix.h:165
static SquareMatrixDouble NewGaussianSobelY(double sigma)
Create Gaussian Sobel Y derivation mask.
Square double matrix class.
bool IsLowerTriangular() const
Check if matrix is lower triangular.
std::multimap< double, MatrixDouble > MakeSpectralEigensystem() const
Perform diagonalization for 2x2 symmetric matrix.
SquareMatrixDouble MakeMinor(size_t r, size_t c) const
Minor matrix.
A character string class.
Definition: CRNStringUTF8.h:49
std::pair< size_t, size_t > Argmax() const
Definition: CRNMatrix.h:451
#define MULT
double Cofactor(size_t r, size_t c) const
Cofactor.
double Determinant() const
Determinant.
SquareMatrixDouble MakeCholesky() const
Get the lower triangular factor in Cholesky decomposition.
double DiagonalProduct() const
Diagonal product.
#define CRN_BEGIN_CLASS_CONSTRUCTOR(classname)
Defines a class constructor.
Definition: CRNObject.h:185
double Pythagoras(double a, double b) noexcept
Computes sqrt(a²+b²) without destructive underflow or overflow.
Definition: CRNMath.h:84
void SwapRows(size_t r1, size_t r2)
Swap two rows in matrix.
Definition: CRNMatrix.h:385