1 // Copyright (C) 2009 Davis E. King (davis@dlib.net) 2 // License: Boost Software License See LICENSE.txt for the full license. 3 #undef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_ 4 #ifdef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_ 5 6 #include "matrix_abstract.h" 7 #include <complex> 8 9 namespace dlib 10 { 11 12 // ---------------------------------------------------------------------------------------- 13 // ---------------------------------------------------------------------------------------- 14 // Global linear algebra functions 15 // ---------------------------------------------------------------------------------------- 16 // ---------------------------------------------------------------------------------------- 17 18 const matrix_exp::matrix_type inv ( 19 const matrix_exp& m 20 ); 21 /*! 22 requires 23 - m is a square matrix 24 ensures 25 - returns the inverse of m 26 (Note that if m is singular or so close to being singular that there 27 is a lot of numerical error then the returned matrix will be bogus. 28 You can check by seeing if m*inv(m) is an identity matrix) 29 !*/ 30 31 // ---------------------------------------------------------------------------------------- 32 33 const matrix pinv ( 34 const matrix_exp& m, 35 double tol = 0 36 ); 37 /*! 38 requires 39 - tol >= 0 40 ensures 41 - returns the Moore-Penrose pseudoinverse of m. 42 - The returned matrix has m.nc() rows and m.nr() columns. 43 - if (tol == 0) then 44 - singular values less than max(m.nr(),m.nc()) times the machine epsilon 45 times the largest singular value are ignored. 46 - else 47 - singular values less than tol*max(singular value in m) are ignored. 48 !*/ 49 50 // ---------------------------------------------------------------------------------------- 51 52 void svd ( 53 const matrix_exp& m, 54 matrix<matrix_exp::type>& u, 55 matrix<matrix_exp::type>& w, 56 matrix<matrix_exp::type>& v 57 ); 58 /*! 59 ensures 60 - computes the singular value decomposition of m 61 - m == #u*#w*trans(#v) 62 - trans(#u)*#u == identity matrix 63 - trans(#v)*#v == identity matrix 64 - diag(#w) == the singular values of the matrix m in no 65 particular order. All non-diagonal elements of #w are 66 set to 0. 67 - #u.nr() == m.nr() 68 - #u.nc() == m.nc() 69 - #w.nr() == m.nc() 70 - #w.nc() == m.nc() 71 - #v.nr() == m.nc() 72 - #v.nc() == m.nc() 73 - if DLIB_USE_LAPACK is #defined then the xGESVD routine 74 from LAPACK is used to compute the SVD. 75 !*/ 76 77 // ---------------------------------------------------------------------------------------- 78 79 long svd2 ( 80 bool withu, 81 bool withv, 82 const matrix_exp& m, 83 matrix<matrix_exp::type>& u, 84 matrix<matrix_exp::type>& w, 85 matrix<matrix_exp::type>& v 86 ); 87 /*! 88 requires 89 - m.nr() >= m.nc() 90 ensures 91 - computes the singular value decomposition of matrix m 92 - m == subm(#u,get_rect(m))*diagm(#w)*trans(#v) 93 - trans(#u)*#u == identity matrix 94 - trans(#v)*#v == identity matrix 95 - #w == the singular values of the matrix m in no 96 particular order. 97 - #u.nr() == m.nr() 98 - #u.nc() == m.nr() 99 - #w.nr() == m.nc() 100 - #w.nc() == 1 101 - #v.nr() == m.nc() 102 - #v.nc() == m.nc() 103 - if (widthu == false) then 104 - ignore the above regarding #u, it isn't computed and its 105 output state is undefined. 106 - if (widthv == false) then 107 - ignore the above regarding #v, it isn't computed and its 108 output state is undefined. 109 - returns an error code of 0, if no errors and 'k' if we fail to 110 converge at the 'kth' singular value. 111 - if (DLIB_USE_LAPACK is #defined) then 112 - if (withu == withv) then 113 - the xGESDD routine from LAPACK is used to compute the SVD. 114 - else 115 - the xGESVD routine from LAPACK is used to compute the SVD. 116 !*/ 117 118 // ---------------------------------------------------------------------------------------- 119 120 void svd3 ( 121 const matrix_exp& m, 122 matrix<matrix_exp::type>& u, 123 matrix<matrix_exp::type>& w, 124 matrix<matrix_exp::type>& v 125 ); 126 /*! 127 ensures 128 - computes the singular value decomposition of m 129 - m == #u*diagm(#w)*trans(#v) 130 - trans(#u)*#u == identity matrix 131 - trans(#v)*#v == identity matrix 132 - #w == the singular values of the matrix m in no 133 particular order. 134 - #u.nr() == m.nr() 135 - #u.nc() == m.nc() 136 - #w.nr() == m.nc() 137 - #w.nc() == 1 138 - #v.nr() == m.nc() 139 - #v.nc() == m.nc() 140 - if DLIB_USE_LAPACK is #defined then the xGESVD routine 141 from LAPACK is used to compute the SVD. 142 !*/ 143 144 // ---------------------------------------------------------------------------------------- 145 146 template < 147 typename T 148 > 149 void svd_fast ( 150 const matrix<T>& A, 151 matrix<T>& u, 152 matrix<T>& w, 153 matrix<T>& v, 154 unsigned long l, 155 unsigned long q = 1 156 ); 157 /*! 158 requires 159 - l > 0 160 - A.size() > 0 161 (i.e. A can't be an empty matrix) 162 ensures 163 - computes the singular value decomposition of A. 164 - Lets define some constants we use to document the behavior of svd_fast(): 165 - Let m = A.nr() 166 - Let n = A.nc() 167 - Let k = min(l, min(m,n)) 168 - Therefore, A represents an m by n matrix and svd_fast() is designed 169 to find a rank-k representation of it. 170 - if (the rank of A is <= k) then 171 - A == #u*diagm(#w)*trans(#v) 172 - else 173 - A is approximated by #u*diagm(#w)*trans(#v) 174 (i.e. In this case A can't be represented with a rank-k matrix, so the 175 matrix you get by trying to reconstruct A from the output of the SVD is 176 not exactly the same.) 177 - trans(#u)*#u == identity matrix 178 - trans(#v)*#v == identity matrix 179 - #w == the top k singular values of the matrix A (in no particular order). 180 - #u.nr() == m 181 - #u.nc() == k 182 - #w.nr() == k 183 - #w.nc() == 1 184 - #v.nr() == n 185 - #v.nc() == k 186 - This function implements the randomized subspace iteration defined in the 187 algorithm 4.4 and 5.1 boxes of the paper: 188 Finding Structure with Randomness: Probabilistic Algorithms for 189 Constructing Approximate Matrix Decompositions by Halko et al. 190 Therefore, it is very fast and suitable for use with very large matrices. 191 Moreover, q is the number of subspace iterations performed. Larger 192 values of q might increase the accuracy of the solution but the default 193 value should be good for many problems. 194 !*/ 195 196 // ---------------------------------------------------------------------------------------- 197 198 template < 199 typename sparse_vector_type, 200 typename T 201 > 202 void svd_fast ( 203 const std::vector<sparse_vector_type>& A, 204 matrix<T>& u, 205 matrix<T>& w, 206 matrix<T>& v, 207 unsigned long l, 208 unsigned long q = 1 209 ); 210 /*! 211 requires 212 - A contains a set of sparse vectors. See dlib/svm/sparse_vector_abstract.h 213 for a definition of what constitutes a sparse vector. 214 - l > 0 215 - max_index_plus_one(A) > 0 216 (i.e. A can't be an empty matrix) 217 ensures 218 - computes the singular value decomposition of A. In this case, we interpret A 219 as a matrix of A.size() rows, where each row is defined by a sparse vector. 220 - Lets define some constants we use to document the behavior of svd_fast(): 221 - Let m = A.size() 222 - Let n = max_index_plus_one(A) 223 - Let k = min(l, min(m,n)) 224 - Therefore, A represents an m by n matrix and svd_fast() is designed 225 to find a rank-k representation of it. 226 - if (the rank of A is <= k) then 227 - A == #u*diagm(#w)*trans(#v) 228 - else 229 - A is approximated by #u*diagm(#w)*trans(#v) 230 (i.e. In this case A can't be represented with a rank-k matrix, so the 231 matrix you get by trying to reconstruct A from the output of the SVD is 232 not exactly the same.) 233 - trans(#u)*#u == identity matrix 234 - trans(#v)*#v == identity matrix 235 - #w == the top k singular values of the matrix A (in no particular order). 236 - #u.nr() == m 237 - #u.nc() == k 238 - #w.nr() == k 239 - #w.nc() == 1 240 - #v.nr() == n 241 - #v.nc() == k 242 - This function implements the randomized subspace iteration defined in the 243 algorithm 4.4 and 5.1 boxes of the paper: 244 Finding Structure with Randomness: Probabilistic Algorithms for 245 Constructing Approximate Matrix Decompositions by Halko et al. 246 Therefore, it is very fast and suitable for use with very large matrices. 247 Moreover, q is the number of subspace iterations performed. Larger 248 values of q might increase the accuracy of the solution but the default 249 value should be good for many problems. 250 !*/ 251 252 template < 253 typename sparse_vector_type, 254 typename T 255 > 256 void svd_fast ( 257 const std::vector<sparse_vector_type>& A, 258 matrix<T>& w, 259 matrix<T>& v, 260 unsigned long l, 261 unsigned long q = 1 262 ); 263 /*! 264 This function is identical to the above svd_fast() except it doesn't compute u. 265 !*/ 266 267 // ---------------------------------------------------------------------------------------- 268 269 template < 270 typename T, 271 long NR, 272 long NC, 273 typename MM, 274 typename L 275 > 276 void orthogonalize ( 277 matrix<T,NR,NC,MM,L>& m 278 ); 279 /*! 280 requires 281 - m.nr() >= m.nc() 282 - m.size() > 0 283 ensures 284 - #m == an orthogonal matrix with the same dimensions as m. In particular, 285 the columns of #m have the same span as the columns of m. 286 - trans(#m)*#m == identity matrix 287 - This function is just shorthand for computing the QR decomposition of m 288 and then storing the Q factor into #m. 289 !*/ 290 291 // ---------------------------------------------------------------------------------------- 292 293 const matrix real_eigenvalues ( 294 const matrix_exp& m 295 ); 296 /*! 297 requires 298 - m.nr() == m.nc() 299 - matrix_exp::type == float or double 300 ensures 301 - returns a matrix E such that: 302 - E.nr() == m.nr() 303 - E.nc() == 1 304 - E contains the real part of all eigenvalues of the matrix m. 305 (note that the eigenvalues are not sorted) 306 !*/ 307 308 // ---------------------------------------------------------------------------------------- 309 310 const matrix_exp::type det ( 311 const matrix_exp& m 312 ); 313 /*! 314 requires 315 - m is a square matrix 316 ensures 317 - returns the determinant of m 318 !*/ 319 320 // ---------------------------------------------------------------------------------------- 321 322 const matrix_exp::type trace ( 323 const matrix_exp& m 324 ); 325 /*! 326 requires 327 - m is a square matrix 328 ensures 329 - returns the trace of m 330 (i.e. returns sum(diag(m))) 331 !*/ 332 333 // ---------------------------------------------------------------------------------------- 334 335 const matrix_exp::matrix_type chol ( 336 const matrix_exp& A 337 ); 338 /*! 339 requires 340 - A is a square matrix 341 ensures 342 - if (A has a Cholesky Decomposition) then 343 - returns the decomposition of A. That is, returns a matrix L 344 such that L*trans(L) == A. L will also be lower triangular. 345 - else 346 - returns a matrix with the same dimensions as A but it 347 will have a bogus value. I.e. it won't be a decomposition. 348 In this case the algorithm returns a partial decomposition. 349 - You can tell when chol fails by looking at the lower right 350 element of the returned matrix. If it is 0 then it means 351 A does not have a cholesky decomposition. 352 353 - If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF 354 is used to compute the cholesky decomposition. 355 !*/ 356 357 // ---------------------------------------------------------------------------------------- 358 359 const matrix_exp::matrix_type inv_lower_triangular ( 360 const matrix_exp& A 361 ); 362 /*! 363 requires 364 - A is a square matrix 365 ensures 366 - if (A is lower triangular) then 367 - returns the inverse of A. 368 - else 369 - returns a matrix with the same dimensions as A but it 370 will have a bogus value. I.e. it won't be an inverse. 371 !*/ 372 373 // ---------------------------------------------------------------------------------------- 374 375 const matrix_exp::matrix_type inv_upper_triangular ( 376 const matrix_exp& A 377 ); 378 /*! 379 requires 380 - A is a square matrix 381 ensures 382 - if (A is upper triangular) then 383 - returns the inverse of A. 384 - else 385 - returns a matrix with the same dimensions as A but it 386 will have a bogus value. I.e. it won't be an inverse. 387 !*/ 388 389 // ---------------------------------------------------------------------------------------- 390 // ---------------------------------------------------------------------------------------- 391 // Matrix decomposition classes 392 // ---------------------------------------------------------------------------------------- 393 // ---------------------------------------------------------------------------------------- 394 395 template < 396 typename matrix_exp_type 397 > 398 class lu_decomposition 399 { 400 /*! 401 REQUIREMENTS ON matrix_exp_type 402 must be some kind of matrix expression as defined in the 403 dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object) 404 The matrix type must also contain float or double values. 405 406 WHAT THIS OBJECT REPRESENTS 407 This object represents something that can compute an LU 408 decomposition of a real valued matrix. That is, for any 409 matrix A it computes matrices L, U, and a pivot vector P such 410 that rowm(A,P) == L*U. 411 412 The LU decomposition with pivoting always exists, even if the matrix is 413 singular, so the constructor will never fail. The primary use of the 414 LU decomposition is in the solution of square systems of simultaneous 415 linear equations. This will fail if is_singular() returns true (or 416 if A is very nearly singular). 417 418 If DLIB_USE_LAPACK is defined then the LAPACK routine xGETRF 419 is used to compute the LU decomposition. 420 !*/ 421 422 public: 423 424 const static long NR = matrix_exp_type::NR; 425 const static long NC = matrix_exp_type::NC; 426 typedef typename matrix_exp_type::type type; 427 typedef typename matrix_exp_type::mem_manager_type mem_manager_type; 428 typedef typename matrix_exp_type::layout_type layout_type; 429 430 typedef matrix<type,0,0,mem_manager_type,layout_type> matrix_type; 431 typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type; 432 typedef matrix<long,NR,1,mem_manager_type,layout_type> pivot_column_vector_type; 433 434 template <typename EXP> 435 lu_decomposition ( 436 const matrix_exp<EXP> &A 437 ); 438 /*! 439 requires 440 - EXP::type == lu_decomposition::type 441 - A.size() > 0 442 ensures 443 - #nr() == A.nr() 444 - #nc() == A.nc() 445 - #is_square() == (A.nr() == A.nc()) 446 - computes the LU factorization of the given A matrix. 447 !*/ 448 449 bool is_square ( 450 ) const; 451 /*! 452 ensures 453 - if (the input A matrix was a square matrix) then 454 - returns true 455 - else 456 - returns false 457 !*/ 458 459 bool is_singular ( 460 ) const; 461 /*! 462 requires 463 - is_square() == true 464 ensures 465 - if (the input A matrix is singular) then 466 - returns true 467 - else 468 - returns false 469 !*/ 470 471 long nr( 472 ) const; 473 /*! 474 ensures 475 - returns the number of rows in the input matrix 476 !*/ 477 478 long nc( 479 ) const; 480 /*! 481 ensures 482 - returns the number of columns in the input matrix 483 !*/ 484 485 const matrix_type get_l ( 486 ) const; 487 /*! 488 ensures 489 - returns the lower triangular L factor of the LU factorization. 490 - L.nr() == nr() 491 - L.nc() == min(nr(),nc()) 492 !*/ 493 494 const matrix_type get_u ( 495 ) const; 496 /*! 497 ensures 498 - returns the upper triangular U factor of the LU factorization. 499 - U.nr() == min(nr(),nc()) 500 - U.nc() == nc() 501 !*/ 502 503 const pivot_column_vector_type& get_pivot ( 504 ) const; 505 /*! 506 ensures 507 - returns the pivot permutation vector. That is, 508 if A is the input matrix then this function 509 returns a vector P such that: 510 - rowm(A,P) == get_l()*get_u() 511 - P.nr() == A.nr() 512 !*/ 513 514 type det ( 515 ) const; 516 /*! 517 requires 518 - is_square() == true 519 ensures 520 - computes and returns the determinant of the input 521 matrix using LU factors. 522 !*/ 523 524 template <typename EXP> 525 const matrix_type solve ( 526 const matrix_exp<EXP> &B 527 ) const; 528 /*! 529 requires 530 - EXP::type == lu_decomposition::type 531 - is_square() == true 532 - B.nr() == nr() 533 ensures 534 - Let A denote the input matrix to this class's constructor. 535 Then this function solves A*X == B for X and returns X. 536 - Note that if A is singular (or very close to singular) then 537 the X returned by this function won't fit A*X == B very well (if at all). 538 !*/ 539 540 }; 541 542 // ---------------------------------------------------------------------------------------- 543 544 template < 545 typename matrix_exp_type 546 > 547 class cholesky_decomposition 548 { 549 /*! 550 REQUIREMENTS ON matrix_exp_type 551 must be some kind of matrix expression as defined in the 552 dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object) 553 The matrix type must also contain float or double values. 554 555 WHAT THIS OBJECT REPRESENTS 556 This object represents something that can compute a cholesky 557 decomposition of a real valued matrix. That is, for any 558 symmetric, positive definite matrix A, it computes a lower 559 triangular matrix L such that A == L*trans(L). 560 561 If the matrix is not symmetric or positive definite, the function 562 computes only a partial decomposition. This can be tested with 563 the is_spd() flag. 564 565 If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF 566 is used to compute the cholesky decomposition. 567 !*/ 568 569 public: 570 571 const static long NR = matrix_exp_type::NR; 572 const static long NC = matrix_exp_type::NC; 573 typedef typename matrix_exp_type::type type; 574 typedef typename matrix_exp_type::mem_manager_type mem_manager_type; 575 typedef typename matrix_exp_type::layout_type layout_type; 576 577 typedef typename matrix_exp_type::matrix_type matrix_type; 578 typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type; 579 580 template <typename EXP> 581 cholesky_decomposition( 582 const matrix_exp<EXP>& A 583 ); 584 /*! 585 requires 586 - EXP::type == cholesky_decomposition::type 587 - A.size() > 0 588 - A.nr() == A.nc() 589 (i.e. A must be a square matrix) 590 ensures 591 - if (A is symmetric positive-definite) then 592 - #is_spd() == true 593 - Constructs a lower triangular matrix L, such that L*trans(L) == A. 594 and #get_l() == L 595 - else 596 - #is_spd() == false 597 !*/ 598 599 bool is_spd( 600 ) const; 601 /*! 602 ensures 603 - if (the input matrix was symmetric positive-definite) then 604 - returns true 605 - else 606 - returns false 607 !*/ 608 609 const matrix_type& get_l( 610 ) const; 611 /*! 612 ensures 613 - returns the lower triangular factor, L, such that L*trans(L) == A 614 (where A is the input matrix to this class's constructor) 615 - Note that if A is not symmetric positive definite or positive semi-definite 616 then the equation L*trans(L) == A won't hold. 617 !*/ 618 619 template <typename EXP> 620 const matrix solve ( 621 const matrix_exp<EXP>& B 622 ) const; 623 /*! 624 requires 625 - EXP::type == cholesky_decomposition::type 626 - B.nr() == get_l().nr() 627 (i.e. the number of rows in B must match the number of rows in the 628 input matrix A) 629 ensures 630 - Let A denote the input matrix to this class's constructor. Then 631 this function solves A*X = B for X and returns X. 632 - Note that if is_spd() == false or A was really close to being 633 non-SPD then the solver will fail to find an accurate solution. 634 !*/ 635 636 }; 637 638 // ---------------------------------------------------------------------------------------- 639 640 template < 641 typename matrix_exp_type 642 > 643 class qr_decomposition 644 { 645 /*! 646 REQUIREMENTS ON matrix_exp_type 647 must be some kind of matrix expression as defined in the 648 dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object) 649 The matrix type must also contain float or double values. 650 651 WHAT THIS OBJECT REPRESENTS 652 This object represents something that can compute a classical 653 QR decomposition of an m-by-n real valued matrix A with m >= n. 654 655 The QR decomposition is an m-by-n orthogonal matrix Q and an 656 n-by-n upper triangular matrix R so that A == Q*R. The QR decomposition 657 always exists, even if the matrix does not have full rank, so the 658 constructor will never fail. The primary use of the QR decomposition 659 is in the least squares solution of non-square systems of simultaneous 660 linear equations. This will fail if is_full_rank() returns false or 661 A is very nearly not full rank. 662 663 The Q and R factors can be retrieved via the get_q() and get_r() 664 methods. Furthermore, a solve() method is provided to find the 665 least squares solution of Ax=b using the QR factors. 666 667 If DLIB_USE_LAPACK is #defined then the xGEQRF routine 668 from LAPACK is used to compute the QR decomposition. 669 !*/ 670 671 public: 672 673 const static long NR = matrix_exp_type::NR; 674 const static long NC = matrix_exp_type::NC; 675 typedef typename matrix_exp_type::type type; 676 typedef typename matrix_exp_type::mem_manager_type mem_manager_type; 677 typedef typename matrix_exp_type::layout_type layout_type; 678 679 typedef matrix<type,0,0,mem_manager_type,layout_type> matrix_type; 680 681 template <typename EXP> 682 qr_decomposition( 683 const matrix_exp<EXP>& A 684 ); 685 /*! 686 requires 687 - EXP::type == qr_decomposition::type 688 - A.nr() >= A.nc() 689 - A.size() > 0 690 ensures 691 - #nr() == A.nr() 692 - #nc() == A.nc() 693 - computes the QR decomposition of the given A matrix. 694 !*/ 695 696 bool is_full_rank( 697 ) const; 698 /*! 699 ensures 700 - if (the input A matrix had full rank) then 701 - returns true 702 - else 703 - returns false 704 !*/ 705 706 long nr( 707 ) const; 708 /*! 709 ensures 710 - returns the number of rows in the input matrix 711 !*/ 712 713 long nc( 714 ) const; 715 /*! 716 ensures 717 - returns the number of columns in the input matrix 718 !*/ 719 720 const matrix_type get_r ( 721 ) const; 722 /*! 723 ensures 724 - returns a matrix R such that: 725 - R is the upper triangular factor, R, of the QR factorization 726 - get_q()*R == input matrix A 727 - R.nr() == nc() 728 - R.nc() == nc() 729 !*/ 730 731 const matrix_type get_q ( 732 ) const; 733 /*! 734 ensures 735 - returns a matrix Q such that: 736 - Q is the economy-sized orthogonal factor Q from the QR 737 factorization. 738 - trans(Q)*Q == identity matrix 739 - Q*get_r() == input matrix A 740 - Q.nr() == nr() 741 - Q.nc() == nc() 742 !*/ 743 744 void get_q ( 745 matrix_type& Q 746 ) const; 747 /*! 748 ensures 749 - #Q == get_q() 750 - This function exists to allow a user to get the Q matrix without the 751 overhead of returning a matrix by value. 752 !*/ 753 754 template <typename EXP> 755 const matrix_type solve ( 756 const matrix_exp<EXP>& B 757 ) const; 758 /*! 759 requires 760 - EXP::type == qr_decomposition::type 761 - B.nr() == nr() 762 ensures 763 - Let A denote the input matrix to this class's constructor. 764 Then this function finds the least squares solution to the equation A*X = B 765 and returns X. X has the following properties: 766 - X is the matrix that minimizes the two norm of A*X-B. That is, it 767 minimizes sum(squared(A*X - B)). 768 - X.nr() == nc() 769 - X.nc() == B.nc() 770 - Note that this function will fail to output a good solution if is_full_rank() == false 771 or the A matrix is close to not being full rank. 772 !*/ 773 774 }; 775 776 // ---------------------------------------------------------------------------------------- 777 778 template < 779 typename matrix_exp_type 780 > 781 class eigenvalue_decomposition 782 { 783 /*! 784 REQUIREMENTS ON matrix_exp_type 785 must be some kind of matrix expression as defined in the 786 dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object) 787 The matrix type must also contain float or double values. 788 789 WHAT THIS OBJECT REPRESENTS 790 This object represents something that can compute an eigenvalue 791 decomposition of a real valued matrix. So it gives 792 you the set of eigenvalues and eigenvectors for a matrix. 793 794 Let A denote the input matrix to this object's constructor. Then 795 what this object does is it finds two matrices, D and V, such that 796 - A*V == V*D 797 Where V is a square matrix that contains all the eigenvectors 798 of the A matrix (each column of V is an eigenvector) and 799 D is a diagonal matrix containing the eigenvalues of A. 800 801 802 It is important to note that if A is symmetric or non-symmetric you 803 get somewhat different results. If A is a symmetric matrix (i.e. A == trans(A)) 804 then: 805 - All the eigenvalues and eigenvectors of A are real numbers. 806 - Because of this there isn't really any point in using the 807 part of this class's interface that returns complex matrices. 808 All you need are the get_real_eigenvalues() and 809 get_pseudo_v() functions. 810 - V*trans(V) should be equal to the identity matrix. That is, all the 811 eigenvectors in V should be orthonormal. 812 - So A == V*D*trans(V) 813 - If DLIB_USE_LAPACK is #defined then this object uses the xSYEVR LAPACK 814 routine. 815 816 On the other hand, if A is not symmetric then: 817 - Some of the eigenvalues and eigenvectors might be complex numbers. 818 - An eigenvalue is complex if and only if its corresponding eigenvector 819 is complex. So you can check for this case by just checking 820 get_imag_eigenvalues() to see if any values are non-zero. You don't 821 have to check the V matrix as well. 822 - V*trans(V) won't be equal to the identity matrix but it is usually 823 invertible. So A == V*D*inv(V) is usually a valid statement but 824 A == V*D*trans(V) won't be. 825 - If DLIB_USE_LAPACK is #defined then this object uses the xGEEV LAPACK 826 routine. 827 !*/ 828 829 public: 830 831 const static long NR = matrix_exp_type::NR; 832 const static long NC = matrix_exp_type::NC; 833 typedef typename matrix_exp_type::type type; 834 typedef typename matrix_exp_type::mem_manager_type mem_manager_type; 835 typedef typename matrix_exp_type::layout_type layout_type; 836 837 typedef typename matrix_exp_type::matrix_type matrix_type; 838 typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type; 839 840 typedef matrix<std::complex<type>,0,0,mem_manager_type,layout_type> complex_matrix_type; 841 typedef matrix<std::complex<type>,NR,1,mem_manager_type,layout_type> complex_column_vector_type; 842 843 844 template <typename EXP> 845 eigenvalue_decomposition( 846 const matrix_exp<EXP>& A 847 ); 848 /*! 849 requires 850 - A.nr() == A.nc() 851 - A.size() > 0 852 - EXP::type == eigenvalue_decomposition::type 853 ensures 854 - #dim() == A.nr() 855 - computes the eigenvalue decomposition of A. 856 - #get_eigenvalues() == the eigenvalues of A 857 - #get_v() == all the eigenvectors of A 858 !*/ 859 860 template <typename EXP> 861 eigenvalue_decomposition( 862 const matrix_op<op_make_symmetric<EXP> >& A 863 ); 864 /*! 865 requires 866 - A.nr() == A.nc() 867 - A.size() > 0 868 - EXP::type == eigenvalue_decomposition::type 869 ensures 870 - #dim() == A.nr() 871 - computes the eigenvalue decomposition of the symmetric matrix A. Does so 872 using a method optimized for symmetric matrices. 873 - #get_eigenvalues() == the eigenvalues of A 874 - #get_v() == all the eigenvectors of A 875 - moreover, since A is symmetric there won't be any imaginary eigenvalues. So 876 we will have: 877 - #get_imag_eigenvalues() == 0 878 - #get_real_eigenvalues() == the eigenvalues of A 879 - #get_pseudo_v() == all the eigenvectors of A 880 - diagm(#get_real_eigenvalues()) == #get_pseudo_d() 881 882 Note that the symmetric matrix operator is created by the 883 dlib::make_symmetric() function. This function simply reflects 884 the lower triangular part of a square matrix into the upper triangular 885 part to create a symmetric matrix. It can also be used to denote that a 886 matrix is already symmetric using the C++ type system. 887 !*/ 888 889 long dim ( 890 ) const; 891 /*! 892 ensures 893 - dim() == the number of rows/columns in the input matrix A 894 !*/ 895 896 const complex_column_vector_type get_eigenvalues ( 897 ) const; 898 /*! 899 ensures 900 - returns diag(get_d()). That is, returns a 901 vector that contains the eigenvalues of the input 902 matrix. 903 - the returned vector has dim() rows 904 - the eigenvalues are not sorted in any particular way 905 !*/ 906 907 const column_vector_type& get_real_eigenvalues ( 908 ) const; 909 /*! 910 ensures 911 - returns the real parts of the eigenvalues. That is, 912 returns real(get_eigenvalues()) 913 - the returned vector has dim() rows 914 - the eigenvalues are not sorted in any particular way 915 !*/ 916 917 const column_vector_type& get_imag_eigenvalues ( 918 ) const; 919 /*! 920 ensures 921 - returns the imaginary parts of the eigenvalues. That is, 922 returns imag(get_eigenvalues()) 923 - the returned vector has dim() rows 924 - the eigenvalues are not sorted in any particular way 925 !*/ 926 927 const complex_matrix_type get_v ( 928 ) const; 929 /*! 930 ensures 931 - returns the eigenvector matrix V that is 932 dim() rows by dim() columns 933 - Each column in V is one of the eigenvectors of the input 934 matrix 935 !*/ 936 937 const complex_matrix_type get_d ( 938 ) const; 939 /*! 940 ensures 941 - returns a matrix D such that: 942 - D.nr() == dim() 943 - D.nc() == dim() 944 - diag(D) == get_eigenvalues() 945 (i.e. the diagonal of D contains all the eigenvalues in the input matrix) 946 - all off diagonal elements of D are set to 0 947 !*/ 948 949 const matrix_type& get_pseudo_v ( 950 ) const; 951 /*! 952 ensures 953 - returns a matrix that is dim() rows by dim() columns 954 - Let A denote the input matrix given to this object's constructor. 955 - if (A has any imaginary eigenvalues) then 956 - returns the pseudo-eigenvector matrix V 957 - The matrix V returned by this function is structured such that: 958 - A*V == V*get_pseudo_d() 959 - else 960 - returns the eigenvector matrix V with A's eigenvectors as 961 the columns of V 962 - A*V == V*diagm(get_real_eigenvalues()) 963 !*/ 964 965 const matrix_type get_pseudo_d ( 966 ) const; 967 /*! 968 ensures 969 - The returned matrix is dim() rows by dim() columns 970 - Computes and returns the block diagonal eigenvalue matrix. 971 If the original matrix A is not symmetric, then the eigenvalue 972 matrix D is block diagonal with the real eigenvalues in 1-by-1 973 blocks and any complex eigenvalues, 974 a + i*b, in 2-by-2 blocks, (a, b; -b, a). That is, if the complex 975 eigenvalues look like 976 977 u + iv . . . . . 978 . u - iv . . . . 979 . . a + ib . . . 980 . . . a - ib . . 981 . . . . x . 982 . . . . . y 983 984 Then D looks like 985 986 u v . . . . 987 -v u . . . . 988 . . a b . . 989 . . -b a . . 990 . . . . x . 991 . . . . . y 992 993 This keeps V (The V you get from get_pseudo_v()) a real matrix in both 994 symmetric and non-symmetric cases, and A*V = V*D. 995 - the eigenvalues are not sorted in any particular way 996 !*/ 997 998 }; 999 1000 // ---------------------------------------------------------------------------------------- 1001 1002 } 1003 1004 #endif // DLIB_MATRIx_LA_FUNCTS_ABSTRACT_ 1005 1006