Line data Source code
1 : /** \file eigenLapack.hpp
2 : * \brief Interfaces to Lapack and BLAS for Eigen-like arrays.
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup gen_math_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
12 : //
13 : // This file is part of mxlib.
14 : //
15 : // mxlib is free software: you can redistribute it and/or modify
16 : // it under the terms of the GNU General Public License as published by
17 : // the Free Software Foundation, either version 3 of the License, or
18 : // (at your option) any later version.
19 : //
20 : // mxlib is distributed in the hope that it will be useful,
21 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : // GNU General Public License for more details.
24 : //
25 : // You should have received a copy of the GNU General Public License
26 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
27 : //***********************************************************************//
28 :
29 : #ifndef math_eigenLapack_hpp
30 : #define math_eigenLapack_hpp
31 :
32 : #pragma GCC system_header
33 : #include <Eigen/Dense>
34 :
35 : #include <cstddef>
36 : #include <cmath>
37 : #include <cstdlib>
38 : #include <limits>
39 :
40 : #include "floatUtils.hpp"
41 : #include "templateBLAS.hpp"
42 : #include "templateLapack.hpp"
43 :
44 : #include "../sys/timeUtils.hpp"
45 :
46 : // #include "vectorUtils.hpp"
47 :
48 : #include "../math/plot/gnuPlot.hpp"
49 :
50 : namespace mx
51 : {
52 : namespace math
53 : {
54 :
55 : /// \cond eigenLapack_test_detail
56 : namespace detail
57 : {
58 :
59 : /// Test-only injection points for deterministic eigensolver failure coverage.
60 : /** Production code must leave both hooks null so the normal allocator and LAPACK wrapper are used.
61 : */
62 : template <typename floatT>
63 : struct eigenLapackTestHooks
64 : {
65 : /// Malloc-compatible allocation function used by workspace-growth tests.
66 : using allocatorT = void *(*)( std::size_t );
67 :
68 : /// Function signature of the LAPACK SYEVR wrapper.
69 : using solverT = MXLAPACK_INT ( * )( char,
70 : char,
71 : char,
72 : MXLAPACK_INT,
73 : floatT *,
74 : MXLAPACK_INT,
75 : floatT,
76 : floatT,
77 : MXLAPACK_INT,
78 : MXLAPACK_INT,
79 : floatT,
80 : MXLAPACK_INT *,
81 : floatT *,
82 : floatT *,
83 : MXLAPACK_INT,
84 : MXLAPACK_INT *,
85 : floatT *,
86 : MXLAPACK_INT,
87 : MXLAPACK_INT *,
88 : MXLAPACK_INT );
89 :
90 : /// Optional malloc-compatible test allocator; null selects `::malloc`.
91 : static inline allocatorT allocator{ nullptr };
92 :
93 : /// Optional test eigensolver; null selects `math::syevr<floatT>`.
94 : static inline solverT solver{ nullptr };
95 : };
96 :
97 : /// Allocate or grow a SYEVR workspace buffer without discarding a usable buffer on failure.
98 : template <typename floatT, typename itemT>
99 113 : bool resizeSyevrBuffer( itemT *&buffer, /**< [in,out] workspace buffer */
100 : MXLAPACK_INT &capacity, /**< [in,out] current item capacity */
101 : MXLAPACK_INT requestedItems /**< [in] required item capacity */ )
102 : {
103 113 : if( capacity >= requestedItems )
104 : {
105 10 : return true;
106 : }
107 :
108 103 : void *rawBuffer{ nullptr };
109 103 : if( eigenLapackTestHooks<floatT>::allocator )
110 : {
111 21 : rawBuffer = eigenLapackTestHooks<floatT>::allocator( requestedItems * sizeof( itemT ) );
112 : }
113 : else
114 : {
115 82 : rawBuffer = ::malloc( requestedItems * sizeof( itemT ) );
116 : }
117 :
118 103 : if( rawBuffer == nullptr )
119 : {
120 7 : return false;
121 : }
122 :
123 96 : if( buffer )
124 : {
125 5 : ::free( buffer );
126 : }
127 :
128 96 : buffer = static_cast<itemT *>( rawBuffer );
129 96 : capacity = requestedItems;
130 96 : return true;
131 : }
132 :
133 : /// Release a SYEVR workspace buffer and reset its ownership state.
134 : template <typename itemT>
135 155 : void releaseSyevrBuffer( itemT *&buffer, /**< [in,out] workspace buffer */
136 : MXLAPACK_INT &capacity /**< [out] reset item capacity */ )
137 : {
138 155 : if( buffer )
139 : {
140 91 : ::free( buffer );
141 : }
142 :
143 155 : buffer = nullptr;
144 155 : capacity = 0;
145 155 : }
146 :
147 : /// Invoke the injected test eigensolver or the production LAPACK wrapper.
148 : template <typename floatT>
149 42 : MXLAPACK_INT callSyevr( char JOBZ, /**< [in] eigenvector request */
150 : char RANGE, /**< [in] eigenvalue selection mode */
151 : char UPLO, /**< [in] populated input triangle */
152 : MXLAPACK_INT N, /**< [in] matrix order */
153 : floatT *A, /**< [in,out] input matrix */
154 : MXLAPACK_INT LDA, /**< [in] input leading dimension */
155 : floatT VL, /**< [in] lower value bound */
156 : floatT VU, /**< [in] upper value bound */
157 : MXLAPACK_INT IL, /**< [in] one-based lower index */
158 : MXLAPACK_INT IU, /**< [in] one-based upper index */
159 : floatT ABSTOL, /**< [in] convergence tolerance */
160 : MXLAPACK_INT *M, /**< [out] selected eigenvalue count */
161 : floatT *W, /**< [out] eigenvalues */
162 : floatT *Z, /**< [out] eigenvectors */
163 : MXLAPACK_INT LDZ, /**< [in] eigenvector leading dimension */
164 : MXLAPACK_INT *ISUPPZ, /**< [out] eigenvector support */
165 : floatT *WORK, /**< [in,out] floating workspace */
166 : MXLAPACK_INT LWORK, /**< [in] floating workspace size */
167 : MXLAPACK_INT *IWORK, /**< [in,out] integer workspace */
168 : MXLAPACK_INT LIWORK /**< [in] integer workspace size */ )
169 : {
170 42 : if( eigenLapackTestHooks<floatT>::solver )
171 : {
172 15 : return eigenLapackTestHooks<floatT>::solver( JOBZ,
173 : RANGE,
174 : UPLO,
175 : N,
176 : A,
177 : LDA,
178 : VL,
179 : VU,
180 : IL,
181 : IU,
182 : ABSTOL,
183 : M,
184 : W,
185 : Z,
186 : LDZ,
187 : ISUPPZ,
188 : WORK,
189 : LWORK,
190 : IWORK,
191 15 : LIWORK );
192 : }
193 :
194 27 : return math::syevr<floatT>( JOBZ,
195 : RANGE,
196 : UPLO,
197 : N,
198 : A,
199 : LDA,
200 : VL,
201 : VU,
202 : IL,
203 : IU,
204 : ABSTOL,
205 : M,
206 : W,
207 : Z,
208 : LDZ,
209 : ISUPPZ,
210 : WORK,
211 : LWORK,
212 : IWORK,
213 27 : LIWORK );
214 : }
215 :
216 : } // namespace detail
217 : /// \endcond
218 :
219 : /// Calculates the lower triangular part of the covariance matrix of ims.
220 : /** Uses cblas_ssyrk. cv is resized to ims.cols() X ims.cols().
221 : * Calculates \f$ cv = A^T*A \f$.
222 : *
223 : *
224 : * \tparam eigenT1 is the eigen matrix/array type of cv.
225 : * \tparam eigenT2 is the eigen matrix/array type of ims
226 : *
227 : * \ingroup eigen_lapack
228 : */
229 : template <typename eigenT1, typename eigenT2>
230 1 : void eigenSYRK( eigenT1 &cv, ///< [out] is the eigen matrix/array where to store the result
231 : const eigenT2 &ims ///< [in] is the eigen matrix/array (images as columns) to
232 : ///< calculate the covariance of
233 : )
234 : {
235 1 : cv.resize( ims.cols(), ims.cols() );
236 :
237 3 : math::syrk<typename eigenT1::Scalar>( /*const enum CBLAS_ORDER Order*/ CblasColMajor,
238 : /*const enum CBLAS_UPLO Uplo*/ CblasLower,
239 : /*const enum CBLAS_TRANSPOSE Trans*/ CblasTrans,
240 1 : /*const MXLAPACK_INT N*/ ims.cols(),
241 1 : /*const MXLAPACK_INT K*/ ims.rows(),
242 1 : /*const float alpha*/ 1.0,
243 : /*const float *A*/ ims.data(),
244 1 : /*const MXLAPACK_INT lda*/ ims.rows(),
245 1 : /*const float beta*/ 0.,
246 : /*float *C*/ cv.data(),
247 1 : /*const MXLAPACK_INT ldc*/ cv.rows() );
248 1 : }
249 :
250 : /// A struct to hold the working memory for eigenSYEVR and maintain it between calls if desired.
251 : /** \todo this should have the working memory for the first exploratory call to ?syevr as well.
252 : */
253 : template <typename floatT>
254 : struct syevrMem
255 : {
256 : /// Triangle selection associated with the cached workspace query.
257 : char UPLO{ 'L' };
258 :
259 : /// Matrix order associated with the cached workspace query.
260 : MXLAPACK_INT n{ 0 };
261 :
262 : /// Eigenvalue-range selection associated with the cached workspace query.
263 : char RANGE{ 'A' };
264 :
265 : /// Number of eigenvalues reported by the most recent workspace query.
266 : MXLAPACK_INT numeig{ 0 };
267 :
268 : /// One-based lower eigenvalue index associated with the cached workspace query.
269 : MXLAPACK_INT IL{ 0 };
270 :
271 : /// One-based upper eigenvalue index associated with the cached workspace query.
272 : MXLAPACK_INT IU{ 0 };
273 :
274 : /// Capacity of the eigenvector-support workspace.
275 : MXLAPACK_INT sizeISuppZ{ 0 };
276 :
277 : /// Capacity of the minimum floating-point query workspace.
278 : MXLAPACK_INT sizeMinWork{ 0 };
279 :
280 : /// Capacity of the optimized floating-point workspace.
281 : MXLAPACK_INT sizeWork{ 0 };
282 :
283 : /// Capacity of the minimum integer query workspace.
284 : MXLAPACK_INT sizeMinIWork{ 0 };
285 :
286 : /// Capacity of the optimized integer workspace.
287 : MXLAPACK_INT sizeIWork{ 0 };
288 :
289 : /// LAPACK eigenvector-support workspace owned by this object.
290 : MXLAPACK_INT *iSuppZ{ nullptr };
291 :
292 : /// Minimum floating-point workspace used for LAPACK size queries.
293 : floatT *minWork{ nullptr };
294 :
295 : /// Optimized floating-point LAPACK workspace owned by this object.
296 : floatT *work{ nullptr };
297 :
298 : /// Minimum integer workspace used for LAPACK size queries.
299 : MXLAPACK_INT *minIWork{ nullptr };
300 :
301 : /// Optimized integer LAPACK workspace owned by this object.
302 : MXLAPACK_INT *iWork{ nullptr };
303 :
304 : /// Calculation-type covariance storage used by higher-level eigensolver helpers.
305 : Eigen::Array<floatT, Eigen::Dynamic, Eigen::Dynamic> cvd;
306 :
307 : /// Calculation-type eigenvector storage used by higher-level eigensolver helpers.
308 : Eigen::Array<floatT, Eigen::Dynamic, Eigen::Dynamic> evecsd;
309 :
310 : /// Calculation-type eigenvalue storage used by higher-level eigensolver helpers.
311 : Eigen::Array<floatT, Eigen::Dynamic, Eigen::Dynamic> evalsd;
312 :
313 : /// Construct an empty reusable SYEVR workspace.
314 28 : syevrMem()
315 28 : {
316 28 : }
317 :
318 : /// Copying an owning SYEVR workspace is prohibited.
319 : syevrMem( const syevrMem &other /**< [in] workspace that cannot be copied */ ) = delete;
320 :
321 : /// Copy assignment of an owning SYEVR workspace is prohibited.
322 : syevrMem &operator=( const syevrMem &other /**< [in] workspace that cannot be copied */ ) = delete;
323 :
324 : /// Release every allocation owned by this workspace.
325 28 : ~syevrMem()
326 : {
327 28 : free();
328 28 : }
329 :
330 : /// Release all workspace allocations and reset the cached LAPACK configuration.
331 31 : void free()
332 : {
333 31 : detail::releaseSyevrBuffer( iSuppZ, sizeISuppZ );
334 31 : detail::releaseSyevrBuffer( minWork, sizeMinWork );
335 31 : detail::releaseSyevrBuffer( work, sizeWork );
336 31 : detail::releaseSyevrBuffer( minIWork, sizeMinIWork );
337 31 : detail::releaseSyevrBuffer( iWork, sizeIWork );
338 :
339 31 : UPLO = 'L';
340 31 : n = 0;
341 31 : RANGE = 'A';
342 31 : numeig = 0;
343 31 : IL = 0;
344 31 : IU = 0;
345 31 : }
346 : };
347 :
348 : /// Calculate select eigenvalues and eigenvectors of an Eigen Array
349 : /** Uses the templateLapack wrapper for syevr.
350 : *
351 : * \tparam arrT is the eigen-like type containing the data
352 : *
353 : * \returns -1 for invalid matrix/range geometry or invalid workspace-query results.
354 : * \returns -1000 on an malloc allocation error.
355 : * \returns the return code from syevr (info) otherwise.
356 : *
357 : * \ingroup eigen_lapack
358 : */
359 : template <typename arrT>
360 39 : MXLAPACK_INT eigenSYEVR( arrT &eigvec, /**< [out] will contain the eigenvectors as columns*/
361 : arrT &eigval, /**< [out] will contain the eigenvalues*/
362 : arrT &X, /**< [in] is a square matrix which is either upper
363 : or lower (default) triangular*/
364 : int ev0 = 0, /**< [in] [opt] is the first desired eigenvalue
365 : (default = 0)*/
366 : int ev1 = -1, /**< [in] [opt] exclusive upper bound of the desired eigenvalue indices.
367 : If -1 all eigenvalues are returned.*/
368 : char UPLO = 'L', /**< [in] [opt] specifies whether X is upper ('U')
369 : or lower ('L') triangular.
370 : Default is ('L').*/
371 : syevrMem<typename arrT::Scalar> *mem = 0 /**< [in] [opt] holds the working memory arrays,
372 : can be re-passed to avoid unnecessary
373 : re-allocations*/
374 : )
375 : {
376 : typedef typename arrT::Scalar calcT;
377 :
378 39 : MXLAPACK_INT numeig{ 0 };
379 : MXLAPACK_INT info;
380 39 : char RANGE = 'A';
381 :
382 39 : MXLAPACK_INT n = X.rows();
383 39 : if( n <= 0 || X.cols() != n )
384 : {
385 2 : return -1;
386 : }
387 :
388 37 : if( ev1 < -1 || ( ev1 != -1 && ( ev0 < 0 || ev0 >= ev1 || ev1 > n ) ) )
389 : {
390 4 : return -1;
391 : }
392 :
393 33 : MXLAPACK_INT localMem = 0;
394 :
395 33 : if( mem == 0 )
396 : {
397 9 : mem = new syevrMem<calcT>;
398 9 : localMem = 1;
399 : }
400 :
401 33 : MXLAPACK_INT IL = 1;
402 33 : MXLAPACK_INT IU = n;
403 33 : if( ev1 != -1 )
404 : {
405 13 : RANGE = 'I';
406 13 : IL = ev0 + 1; // This is FORTRAN, after all
407 13 : IU = ev1;
408 : }
409 :
410 33 : eigvec.resize( n, IU - IL + 1 );
411 33 : eigval.resize( n, 1 );
412 :
413 33 : if( UPLO != mem->UPLO || n != mem->n || RANGE != mem->RANGE || mem->IL != IL || mem->IU != IU )
414 : {
415 29 : if( !detail::resizeSyevrBuffer<calcT>( mem->iSuppZ, mem->sizeISuppZ, 2 * n ) ||
416 55 : !detail::resizeSyevrBuffer<calcT>( mem->minWork, mem->sizeMinWork, 26 * n ) ||
417 26 : !detail::resizeSyevrBuffer<calcT>( mem->minIWork, mem->sizeMinIWork, 10 * n ) )
418 : {
419 8 : internal::mxlib_error_report( error_t::allocerr, "malloc failed in eigenSYEVR." );
420 4 : if( localMem )
421 : {
422 1 : delete mem;
423 : }
424 4 : return -1000;
425 : }
426 :
427 : // Query for optimum sizes for workspace
428 25 : info = detail::callSyevr<calcT>( 'V',
429 : RANGE,
430 : UPLO,
431 : n,
432 : X.data(),
433 : n,
434 : 0,
435 : 0,
436 : IL,
437 : IU,
438 : math::lamch<calcT>( 'S' ),
439 : &numeig,
440 : eigval.data(),
441 : eigvec.data(),
442 : n,
443 : mem->iSuppZ,
444 : mem->minWork,
445 : -1,
446 : mem->minIWork,
447 : -1 );
448 :
449 25 : if( info != 0 )
450 : {
451 4 : internal::mxlib_error_report( error_t::lapackerr, "error from SYEVR" );
452 2 : if( localMem )
453 : {
454 1 : delete mem;
455 : }
456 :
457 2 : return info;
458 : }
459 :
460 23 : const calcT queriedWork = mem->minWork[0];
461 23 : const MXLAPACK_INT requestedIWork = mem->minIWork[0];
462 44 : if( !math::isFinite( queriedWork ) || queriedWork < 1 ||
463 19 : static_cast<long double>( queriedWork ) >
464 44 : static_cast<long double>( std::numeric_limits<MXLAPACK_INT>::max() ) ||
465 : requestedIWork < 1 )
466 : {
467 14 : internal::mxlib_error_report( error_t::lapackerr, "invalid workspace sizes returned by SYEVR query" );
468 7 : if( localMem )
469 : {
470 5 : delete mem;
471 : }
472 7 : return -1;
473 : }
474 :
475 16 : const MXLAPACK_INT requestedWork = static_cast<MXLAPACK_INT>( queriedWork );
476 31 : if( !detail::resizeSyevrBuffer<calcT>( mem->work, mem->sizeWork, requestedWork ) ||
477 15 : !detail::resizeSyevrBuffer<calcT>( mem->iWork, mem->sizeIWork, requestedIWork ) )
478 : {
479 6 : internal::mxlib_error_report( error_t::allocerr, "malloc failed in eigenSYEVR." );
480 3 : if( localMem )
481 : {
482 1 : delete mem;
483 : }
484 3 : return -1000;
485 : }
486 :
487 13 : mem->UPLO = UPLO;
488 13 : mem->n = n;
489 13 : mem->RANGE = RANGE;
490 13 : mem->numeig = numeig;
491 13 : mem->IL = IL;
492 13 : mem->IU = IU;
493 : }
494 :
495 : // Now actually do the calculation
496 17 : info = detail::callSyevr<calcT>( 'V',
497 : RANGE,
498 : UPLO,
499 : n,
500 : X.data(),
501 : n,
502 : 0,
503 : 0,
504 : IL,
505 : IU,
506 : math::lamch<calcT>( 'S' ),
507 : &numeig,
508 : eigval.data(),
509 : eigvec.data(),
510 : n,
511 : mem->iSuppZ,
512 : mem->work,
513 : mem->sizeWork,
514 : mem->iWork,
515 : mem->sizeIWork );
516 :
517 : /* Cleanup and exit */
518 :
519 17 : if( localMem )
520 : {
521 1 : delete mem;
522 : }
523 :
524 17 : return info;
525 : }
526 :
527 : /// Calculate the eigenvectors and eigenvalues given a triangular matrix
528 : /** Eigen-decomposition of the matrix is performed using \ref eigenSYEVR().
529 : *
530 : * \tparam evCalcT is the type in which to perform eigen-decomposition, which may be different
531 : * from the input array and output arrays (which must be the same).
532 : * \tparam eigenT is a 2D Eigen-like type
533 : *
534 : * \ingroup eigen_lapack
535 : */
536 : template <typename _evCalcT = double, typename eigenT>
537 14 : MXLAPACK_INT calcEigenVecs( eigenT &evecs, /**< [out] on exit contains the eigen vectors*/
538 : eigenT &evals, /**< [out] on exit contains the eigen vectors*/
539 : eigenT &cv, /**< [in] a lower-triangle (in the Lapack sense) square
540 : covariance matrix.*/
541 : int nVecs = 0, /**< [in] [opt] The maximum number of modes to solve
542 : for. If 0 all modes are solved for.*/
543 : bool normalize = false, /**< [in] [opt] flag specifying whether or not to
544 : normalize the eigenvectors.*/
545 : bool check = false, /**< [in] [opt] flag specifying whether or not to
546 : check the eigenvalues/vectors for
547 : validity. Requires normalize=true.*/
548 : syevrMem<_evCalcT> *mem = 0, /**< [in] [opt] A memory structure which can be
549 : re-used by SYEVR for efficiency.*/
550 : double *t_eigenv = nullptr /**< [out] [opt] if not null, will be filled in
551 : with the time taken to calculate
552 : eigenvalues.*/ )
553 : {
554 : typedef _evCalcT evCalcT;
555 : typedef typename eigenT::Scalar realT;
556 :
557 14 : bool localMem = false;
558 :
559 14 : if( mem == 0 )
560 : {
561 4 : mem = new syevrMem<evCalcT>;
562 4 : localMem = true;
563 : }
564 :
565 14 : if( cv.rows() <= 0 || cv.rows() != cv.cols() )
566 : {
567 3 : std::cerr << "calcEigenVecs: covariance matrix must be non-empty and square\n";
568 3 : if( localMem )
569 : {
570 1 : delete mem;
571 : }
572 3 : return -1;
573 : }
574 :
575 11 : MXLAPACK_INT tNims = cv.rows();
576 :
577 11 : if( nVecs <= 0 || nVecs > tNims )
578 : {
579 6 : nVecs = tNims;
580 : }
581 :
582 11 : if( t_eigenv )
583 : {
584 2 : *t_eigenv = sys::get_curr_time();
585 : }
586 :
587 11 : mem->cvd = cv.template cast<evCalcT>();
588 :
589 : // Calculate eigenvectors and eigenvalues
590 : /* SYEVR sorts eigenvalues in ascending order, so we specifiy the top n_modes
591 : */
592 11 : MXLAPACK_INT info = eigenSYEVR( mem->evecsd, mem->evalsd, mem->cvd, tNims - nVecs, tNims, 'L', mem );
593 :
594 11 : if( info != 0 )
595 : {
596 2 : std::cerr << "calcEigenVecs: eigenSYEVR returned an error (info = " << info << ")\n";
597 2 : if( localMem )
598 : {
599 1 : delete mem;
600 : }
601 :
602 2 : return -1;
603 : }
604 :
605 9 : evecs = mem->evecsd.template cast<realT>();
606 9 : evals = mem->evalsd.template cast<realT>();
607 :
608 9 : if( normalize )
609 : {
610 : // Normalize the eigenvectors
611 5 : if( !check )
612 : {
613 3 : for( MXLAPACK_INT i = 0; i < nVecs; ++i )
614 : {
615 2 : evecs.col( i ) = evecs.col( i ) / sqrt( evals( i ) );
616 : }
617 : }
618 : else // here we check for invalid results and 0 things out
619 : {
620 13 : for( MXLAPACK_INT i = 0; i < nVecs; ++i )
621 : {
622 9 : if( !math::isFinite( evals( i ) ) )
623 : {
624 1 : std::cerr << "got non-finite eigenvalue (# " << i << ")\n";
625 1 : evecs.col( i ).setZero();
626 : }
627 8 : else if( evals( i ) == 0 )
628 : {
629 1 : std::cerr << "got 0 eigenvalue (# " << i << ")\n";
630 1 : evecs.col( i ).setZero();
631 : }
632 7 : else if( evals( i ) < 0 )
633 : {
634 1 : std::cerr << "got < 0 eigenvalue (# " << i << ")\n";
635 1 : evecs.col( i ).setZero();
636 : }
637 : else
638 : {
639 6 : evecs.col( i ) = evecs.col( i ) / sqrt( evals( i ) );
640 : }
641 :
642 37 : for( int r = 0; r < evecs.rows(); ++r )
643 : {
644 29 : if( !math::isFinite( evecs.col( i )( r ) ) )
645 : {
646 1 : std::cerr << "got non-finite eigenvector entry (# " << i << "," << r << ")\n";
647 1 : evecs.col( i ).setZero();
648 1 : break;
649 : }
650 : }
651 : }
652 : }
653 : }
654 :
655 9 : if( t_eigenv )
656 : {
657 2 : *t_eigenv = sys::get_curr_time() - *t_eigenv;
658 : }
659 :
660 9 : if( localMem )
661 : {
662 2 : delete mem;
663 : }
664 :
665 9 : return 0;
666 :
667 : } // calcKLModes
668 :
669 : /// Calculate the K-L modes, or principle components, given a covariance matrix.
670 : /** Eigen-decomposition of the covariance matrix is performed using \ref eigenSYEVR().
671 : *
672 : * \tparam evCalcT is the type in which to perform eigen-decomposition.
673 : * \tparam eigenT is a 2D Eigen-like type
674 : * \tparam eigenT1 is a 2D Eigen-like type.
675 : *
676 : * \ingroup eigen_lapack
677 : */
678 : template <typename _evCalcT = double, typename eigenT, typename eigenT1>
679 5 : MXLAPACK_INT calcKLModes( eigenT &klModes, /**< [out] on exit contains the K-L modes (or P.C.s) */
680 : eigenT &cv, /**< [in] a lower-triangle (in the Lapack sense) square
681 : covariance matrix.*/
682 : const eigenT1 &Rims, /**< [in] The reference data. cv.rows() == Rims.cols().*/
683 : int n_modes = 0, /**< [in] [opt] Tbe maximum number of modes to solve for.
684 : If 0 all modes are solved for.*/
685 : syevrMem<_evCalcT> *mem = 0, /**< [in] [opt] A memory structure which can be re-used
686 : by SYEVR for efficiency.*/
687 : double *t_eigenv = nullptr, /**< [out] [opt] if not null, will be filled in with the time
688 : taken to calculate eigenvalues.*/
689 : double *t_klim = nullptr /**< [out] [opt] if not null, will be filled in with the time
690 : taken to calculate the KL modes.*/)
691 : {
692 : typedef _evCalcT evCalcT;
693 : typedef typename eigenT::Scalar realT;
694 :
695 5 : bool localMem = false;
696 :
697 5 : if( mem == 0 )
698 : {
699 3 : mem = new syevrMem<evCalcT>;
700 3 : localMem = true;
701 : }
702 :
703 5 : if( cv.rows() != Rims.cols() )
704 : {
705 1 : std::cerr << "Covariance matrix - reference image size mismatch in calcKLModes\n";
706 1 : if( localMem )
707 : {
708 1 : delete mem;
709 : }
710 1 : return -1;
711 : }
712 :
713 4 : eigenT evecs, evals;
714 :
715 4 : MXLAPACK_INT tNims = cv.rows();
716 4 : MXLAPACK_INT tNpix = Rims.rows();
717 :
718 4 : if( n_modes <= 0 || n_modes > tNims )
719 : {
720 2 : n_modes = tNims;
721 : }
722 :
723 : // Eigen::Array<evCalcT, Eigen::Dynamic, Eigen::Dynamic> evecsd, evalsd;
724 4 : mem->cvd = cv.template cast<evCalcT>();
725 4 : MXLAPACK_INT info = calcEigenVecs( mem->evecsd, mem->evalsd, mem->cvd, n_modes, true, true, mem, t_eigenv );
726 :
727 4 : if( info != 0 )
728 : {
729 1 : std::cerr << "calckKLModes: eigenSYEVR returned an error (info = " << info << ")\n";
730 1 : if( localMem )
731 : {
732 1 : delete mem;
733 : }
734 1 : return -1;
735 : }
736 :
737 3 : evecs = mem->evecsd.template cast<realT>();
738 3 : evals = mem->evalsd.template cast<realT>();
739 :
740 3 : klModes.resize( n_modes, tNpix );
741 :
742 3 : if( t_klim )
743 : {
744 1 : *t_klim = sys::get_curr_time();
745 : }
746 :
747 : // Now calculate KL images
748 : /*
749 : * KL = E^T * R ==> C = A^T * B
750 : */
751 3 : constexpr realT alpha{ 1 };
752 3 : constexpr realT beta{ 0 };
753 6 : gemm<realT>( CblasColMajor,
754 : CblasTrans,
755 : CblasTrans,
756 : n_modes,
757 : tNpix,
758 : tNims,
759 : alpha,
760 3 : evecs.data(),
761 3 : cv.rows(),
762 : Rims.data(),
763 3 : Rims.rows(),
764 : beta,
765 : klModes.data(),
766 3 : klModes.rows() );
767 :
768 3 : if( t_klim )
769 : {
770 1 : *t_klim = sys::get_curr_time() - *t_klim;
771 : }
772 :
773 3 : if( localMem )
774 : {
775 1 : delete mem;
776 : }
777 :
778 3 : return 0;
779 :
780 4 : } // calcKLModes
781 :
782 : /// Compute the SVD of an Eigen::Array using LAPACK's xgesdd
783 : /** Computes the SVD of A, \f$ A = U S V^T \f$.
784 : *
785 : * \returns 0 on success
786 : * \returns -i on error in ith parameter (from LAPACK xgesdd)
787 : * \returns >0 did not converge (from LAPACK xgesdd)
788 : *
789 : * \tparam dataT is either float or double.
790 : *
791 : * \ingroup eigen_lapack
792 : */
793 : template <typename dataT>
794 0 : MXLAPACK_INT eigenGESDD( Eigen::Array<dataT, -1, -1> &U, ///< [out] the A.rows() x A.rows() left matrix
795 : Eigen::Array<dataT, -1, -1> &S, ///< [out] the A.cols() x 1 matrix of singular values
796 : Eigen::Array<dataT, -1, -1> &VT, /**< [out] the A.cols() x A.cols() right matrix, note this
797 : is the transpose. */
798 : Eigen::Array<dataT, -1, -1> &A ///< [in] the input matrix to be decomposed
799 : )
800 : {
801 0 : char JOBZ = 'A';
802 0 : MXLAPACK_INT M = A.rows();
803 0 : MXLAPACK_INT N = A.cols();
804 0 : MXLAPACK_INT LDA = M;
805 0 : S.resize( N, 1 );
806 0 : U.resize( M, M );
807 0 : MXLAPACK_INT LDU = M;
808 0 : VT.resize( N, N );
809 0 : MXLAPACK_INT LDVT = N;
810 :
811 : dataT wkOpt;
812 0 : MXLAPACK_INT LWORK = -1;
813 :
814 0 : MXLAPACK_INT *IWORK = new MXLAPACK_INT[8 * M];
815 : MXLAPACK_INT INFO;
816 :
817 : INFO =
818 0 : math::gesdd<dataT>( JOBZ, M, N, A.data(), LDA, S.data(), U.data(), LDU, VT.data(), LDVT, &wkOpt, LWORK, IWORK );
819 :
820 0 : LWORK = wkOpt;
821 : // delete WORK;
822 0 : dataT *WORK = new dataT[LWORK];
823 :
824 : INFO =
825 0 : math::gesdd<dataT>( JOBZ, M, N, A.data(), LDA, S.data(), U.data(), LDU, VT.data(), LDVT, WORK, LWORK, IWORK );
826 :
827 0 : delete[] WORK;
828 0 : delete[] IWORK;
829 :
830 0 : return INFO;
831 : }
832 :
833 : #define MX_PINV_NO_INTERACT 0
834 : #define MX_PINV_PLOT 1
835 : #define MX_PINV_ASK 2
836 : #define MX_PINV_ASK_NMODES 4
837 :
838 : /// Calculate the pseudo-inverse of a patrix given its SVD
839 : /** Given the SVD of A, \f$ A = U S V^T \f$, as calculated by eigenGESDD the psuedo-inverse is
840 : * calculated as \f$ A^+ = V S^+ U^T\f$.
841 : *
842 : * The parameter \p interact is intepreted as a bitmask. The values can be
843 : * - `MX_PINV_PLOT`, which will cause a plot to be displayed of the singular values
844 : * - `MX_PINV_ASK`, which will ask the user for a max. condition number using stdin
845 : * - `MX_PINV_ASK_NMODES`, which will ask the user for a max number of modes to include using
846 : * stdin. Overrides MX_PINV_ASK. If \p interact is 0 then no interaction is used and \p
847 : * maxCondition controls the inversion.
848 : *
849 : * \tparam dataT is either float or double.
850 : *
851 : * \ingroup eigen_lapack
852 : */
853 : template <typename dataT>
854 0 : int eigenPseudoInverse( Eigen::Array<dataT, -1, -1> &PInv, ///< [out] The pseudo-inverse of A
855 : dataT &condition, ///< [out] The final condition number.
856 : int &nRejected, ///< [out] The number of eigenvectors rejected
857 : Eigen::Array<dataT, -1, -1> &U, ///< [in] the A.rows() x A.rows() left matrix
858 : Eigen::Array<dataT, -1, -1> &S, ///< [in] the A.cols() x 1 matrix of singular values
859 : Eigen::Array<dataT, -1, -1> &VT, /**< [in] the A.cols() x A.cols() right matrix, note this
860 : is the transpose. */
861 : int minMN, ///< [in] The minimum size of the matrix to invert.
862 : dataT &maxCondition, /**< [in] If \> 0, the maximum condition number desired.
863 : If \<0 the number of modes to keep. Used to
864 : threshold the singular values. Set to 0 to
865 : include all eigenvalues/vectors. Ignored if
866 : interactive. */
867 : dataT alpha = 0, /**< [in] [opt] the Tikhonov regularization value, as
868 : a fraction of the highest singular value.
869 : If alpha < 0, then it is treated as a (positive)
870 : floor (as a fraction of highest singular value)
871 : for the singular values, which is not the same as
872 : Tikhonov (alpha > 0).*/
873 : int interact = MX_PINV_NO_INTERACT /**< [in] [opt] a bitmask controlling interaction.
874 : See above.*/
875 : )
876 : {
877 :
878 0 : dataT Smax = S.maxCoeff();
879 :
880 0 : if( alpha > 0 )
881 : {
882 0 : for( MXLAPACK_INT i = 0; i < S.rows(); ++i )
883 : {
884 0 : S( i ) = ( pow( S( i ), 2 ) + pow( alpha * Smax, 2 ) ) / S( i );
885 : }
886 : }
887 :
888 0 : if( alpha < 0 )
889 : {
890 0 : for( MXLAPACK_INT i = 0; i < S.rows(); ++i )
891 : {
892 0 : S( i ) = S( i ) + -alpha * Smax;
893 : }
894 : }
895 :
896 0 : int modesToReject = 0;
897 0 : if( maxCondition < 0 ) // Rejecting mode numbers
898 : {
899 0 : modesToReject = -maxCondition;
900 :
901 0 : if( modesToReject - 1 < S.rows() )
902 : {
903 0 : maxCondition = Smax / S( modesToReject - 1, 0 );
904 : }
905 : }
906 :
907 0 : if( interact & MX_PINV_PLOT )
908 : {
909 0 : gnuPlot gp;
910 0 : gp.command( "set title \"SVD Singular Values\"" );
911 0 : gp.logy();
912 0 : gp.plot( S.data(), S.rows(), " w lp", "singular values" );
913 0 : }
914 :
915 0 : if( interact & MX_PINV_ASK && !( interact & MX_PINV_ASK_NMODES ) )
916 : {
917 : dataT mine;
918 0 : std::cout << "Maximum singular value: " << Smax << "\n";
919 0 : std::cout << "Minimum singular value: " << S.minCoeff() << "\n";
920 0 : std::cout << "Enter singular value threshold: ";
921 0 : std::cin >> mine;
922 :
923 0 : if( mine > 0 )
924 : {
925 0 : maxCondition = Smax / mine;
926 : }
927 : else
928 : {
929 0 : maxCondition = Smax / S( S.rows() - 1, 0 );
930 : }
931 0 : }
932 0 : else if( interact & MX_PINV_ASK_NMODES )
933 : {
934 : unsigned mine;
935 0 : std::cout << "Maximum singular value: " << Smax << "\n";
936 0 : std::cout << "Minimum singular value: " << S.minCoeff() << "\n";
937 0 : std::cout << "Enter number of modes to keep: ";
938 0 : std::cin >> mine;
939 0 : modesToReject = S.rows() - mine;
940 :
941 0 : if( modesToReject <= 0 || modesToReject > S.rows() )
942 : {
943 0 : modesToReject = 0;
944 : }
945 :
946 0 : maxCondition = -modesToReject;
947 : }
948 :
949 0 : Eigen::Array<dataT, -1, -1> sigma;
950 0 : sigma.resize( S.rows(), S.rows() );
951 0 : sigma.setZero();
952 :
953 0 : nRejected = 0;
954 :
955 0 : if( maxCondition > 0 )
956 : {
957 0 : dataT threshold = 0;
958 :
959 0 : if( maxCondition > 0 )
960 : {
961 0 : threshold = Smax / maxCondition;
962 :
963 0 : condition = 1;
964 :
965 0 : for( MXLAPACK_INT i = 0; i < S.rows(); ++i )
966 : {
967 0 : if( S( i ) >= threshold )
968 : {
969 0 : sigma( i, i ) = 1. / S( i );
970 0 : if( Smax / S( i ) > condition )
971 : {
972 0 : condition = Smax / S( i );
973 : }
974 : }
975 : else
976 : {
977 0 : sigma( i, i ) = 0;
978 0 : ++nRejected;
979 : }
980 : }
981 : }
982 : }
983 : else // rejecting modes
984 : {
985 0 : std::cerr << "rejecting based on modes\n";
986 0 : std::cerr << " modes to reject: " << modesToReject << "\n";
987 0 : for( MXLAPACK_INT i = 0; i < S.rows(); ++i )
988 : {
989 0 : if( i < S.rows() - modesToReject )
990 : {
991 0 : sigma( i, i ) = 1. / S( i );
992 0 : if( Smax / S( i ) > condition )
993 0 : condition = Smax / S( i );
994 : }
995 : else
996 : {
997 0 : sigma( i, i ) = 0;
998 0 : ++nRejected;
999 : }
1000 : }
1001 : }
1002 :
1003 0 : if( interact & MX_PINV_PLOT )
1004 : {
1005 0 : std::vector<dataT> vsig( sigma.rows() );
1006 0 : for( int rr = 0; rr < sigma.rows(); ++rr )
1007 : {
1008 0 : vsig[rr] = sigma( rr, rr );
1009 : }
1010 :
1011 0 : gnuPlot gp;
1012 0 : gp.command( "set title \"Inverted Singular Values\"" );
1013 0 : gp.logy();
1014 0 : gp.plot( vsig.data(), vsig.size(), " w lp", "inverted singular values" );
1015 0 : }
1016 :
1017 0 : if( interact & MX_PINV_ASK || interact & MX_PINV_ASK_NMODES )
1018 : {
1019 : dataT mine;
1020 0 : std::cout << "Modes Rejected: " << nRejected << "\n";
1021 0 : std::cout << "Condition Number: " << condition << "\n";
1022 : }
1023 :
1024 0 : PInv = ( VT.matrix().transpose() * sigma.matrix().transpose() ) *
1025 0 : U.block( 0, 0, U.rows(), minMN ).matrix().transpose();
1026 :
1027 0 : return 0;
1028 0 : }
1029 :
1030 : /// Calculate the pseudo-inverse of a patrix using the SVD
1031 : /** First computes the SVD of A, \f$ A = U S V^T \f$, using eigenGESDD. Then the psuedo-inverse is
1032 : * calculated as \f$ A^+ = V S^+ U^T\f$.
1033 : *
1034 : * The parameter \p interact is intepreted as a bitmask. The values can be
1035 : * - `MX_PINV_PLOT`, which will cause a plot to be displayed of the singular values
1036 : * - `MX_PINV_ASK`, which will ask the user for a max. condition number using stdin
1037 : * - `MX_PINV_ASK_NMODES`, which will ask the user for a max number of modes to include using
1038 : * stdin. Overrides MX_PINV_ASK. If \p interact is 0 then no interaction is used and \p
1039 : * maxCondition controls the inversion.
1040 : *
1041 : * \tparam dataT is either float or double.
1042 : *
1043 : * \ingroup eigen_lapack
1044 : */
1045 : template <typename dataT>
1046 0 : int eigenPseudoInverse( Eigen::Array<dataT, -1, -1> &PInv, ///< [out] The pseudo-inverse of A
1047 : dataT &condition, ///< [out] The final condition number.
1048 : int &nRejected, ///< [out] The number of eigenvectors rejected
1049 : Eigen::Array<dataT, -1, -1> &U, ///< [out] the A.rows() x A.rows() left matrix
1050 : Eigen::Array<dataT, -1, -1> &S, ///< [out] the A.cols() x 1 matrix of singular values
1051 : Eigen::Array<dataT, -1, -1> &VT, /**< [out] the A.cols() x A.cols() right matrix, note this
1052 : is the transpose. */
1053 : Eigen::Array<dataT, -1, -1> &A, ///< [in] The matrix to invert. This will be modified!
1054 : dataT &maxCondition, /**< [in] If \> 0, the maximum condition number desired.
1055 : If \<0 the number of modes to keep. Used to
1056 : threshold the singular values. Set to 0 to
1057 : include all eigenvalues/vectors. Ignored if
1058 : interactive. */
1059 : dataT alpha = 0, /**< [in] [opt] the Tikhonov regularization value, as
1060 : a fraction of the highest singular value.
1061 : If alpha < 0, then it is treated as a (positive)
1062 : floor (as a fraction of highest singular value)
1063 : for the singular values, which is not the same as
1064 : Tikhonov (alpha > 0).*/
1065 : int interact = MX_PINV_NO_INTERACT /**< [in] [opt] a bitmask controlling interaction.
1066 : See above.*/
1067 : )
1068 : {
1069 :
1070 0 : int minMN = std::min( A.rows(), A.cols() );
1071 :
1072 : MXLAPACK_INT info;
1073 0 : info = eigenGESDD( U, S, VT, A );
1074 :
1075 0 : if( info != 0 )
1076 : {
1077 0 : std::cerr << "eigenPseudoInverse: eigenGESDD failed with info = " << info << "\n";
1078 0 : return info;
1079 : }
1080 :
1081 0 : return eigenPseudoInverse( PInv, condition, nRejected, U, S, VT, minMN, maxCondition, alpha, interact );
1082 : }
1083 :
1084 : /// Calculate the pseudo-inverse of a matrix using the SVD
1085 : /** First computes the SVD of A, \f$ A = U S V^T \f$, using eigenGESDD. Then the psuedo-inverse is
1086 : * calculated as \f$ A^+ = V S^+ U^T\f$. This interface does not provide access to U, S and VT.
1087 : *
1088 : * The parameter \p interact is intepreted as a bitmask. The values can be
1089 : * - `MX_PINV_PLOT`, which will cause a plot to be displayed of the singular values
1090 : * - `MX_PINV_ASK`, which will ask the user for a max. condition number using stdin
1091 : * - `MX_PINV_ASK_NMODES`, which will ask the user for a max number of modes to include using
1092 : * stdin. Overrides MX_PINV_ASK. If \p interact is 0 then no interaction is used and maxCondition
1093 : * controls the inversion.
1094 : * *
1095 : * \tparam dataT is either float or double.
1096 : *
1097 : * \overload
1098 : *
1099 : * \ingroup eigen_lapack
1100 : */
1101 : template <typename dataT>
1102 0 : int eigenPseudoInverse( Eigen::Array<dataT, -1, -1> &PInv, ///< [out] The pseudo-inverse of A
1103 : dataT &condition, ///< [out] The final condition number.
1104 : int &nRejected, /**< [out] The number of eigenvectors
1105 : rejected*/
1106 : Eigen::Array<dataT, -1, -1> &A, /**< [in] The matrix to invert, will be
1107 : altered!*/
1108 : dataT &maxCondition, /**< [in] If \> 0, the maximum condition number desired.
1109 : If \<0 the number of modes to keep. Used to
1110 : threshold the singular values. Set to 0 to
1111 : include all eigenvalues/vectors. Ignored if
1112 : interactive.*/
1113 : dataT alpha = 0, /**< [in] [opt] the Tikhonov regularization value,
1114 : as a fraction of the highest singular value. If
1115 : alpha < 0, then it is treated as a (positive)
1116 : floor (as a fraction of highest singular value)
1117 : for the singular values, which is not the same
1118 : as Tikhonov (alpha > 0).*/
1119 : int interact = MX_PINV_NO_INTERACT /**< [in] [opt] a bitmask controlling interaction.
1120 : See above.*/
1121 : )
1122 : {
1123 0 : Eigen::Array<dataT, -1, -1> S, U, VT;
1124 :
1125 0 : return eigenPseudoInverse( PInv, condition, nRejected, U, S, VT, A, maxCondition, alpha, interact );
1126 0 : }
1127 :
1128 : } // namespace math
1129 : } // namespace mx
1130 :
1131 : #endif // math_eigenLapack_hpp
|
