Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TDecompSVD.cxx
Go to the documentation of this file.
1// @(#)root/matrix:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Dec 2003
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TDecompSVD
13 \ingroup Matrix
14
15 Single Value Decomposition class
16
17 For an (m x n) matrix A with m >= n, the singular value decomposition
18 is
19 an (m x m) orthogonal matrix fU, an (m x n) diagonal matrix fS, and
20 an (n x n) orthogonal matrix fV so that A = U*S*V'.
21
22 If the row/column index of A starts at (rowLwb,colLwb) then the
23 decomposed matrices/vectors start at :
24~~~
25 fU : (rowLwb,colLwb)
26 fV : (colLwb,colLwb)
27 fSig : (colLwb)
28~~~
29
30 The diagonal matrix fS is stored in the singular values vector fSig .
31 The singular values, fSig[k] = S[k][k], are ordered so that
32 fSig[0] >= fSig[1] >= ... >= fSig[n-1].
33
34 The singular value decomposition always exists, so the decomposition
35 will (as long as m >=n) never fail. If m < n, the user should add
36 sufficient zero rows to A , so that m == n
37
38 Here fTol is used to set the threshold on the minimum allowed value
39 of the singular values:
40 min_singular = fTol*max(fSig[i])
41*/
42
43#include "TDecompSVD.h"
44#include "TMath.h"
45#include "TArrayD.h"
46
47
48////////////////////////////////////////////////////////////////////////////////
49/// Constructor for (nrows x ncols) matrix
50
52{
53 if (nrows < ncols) {
54 Error("TDecompSVD(Int_t,Int_t","matrix rows should be >= columns");
55 return;
56 }
59 fV.ResizeTo(nrows,ncols); // In the end we only need the nColxnCol part
60}
61
62////////////////////////////////////////////////////////////////////////////////
63/// Constructor for ([row_lwb..row_upb] x [col_lwb..col_upb]) matrix
64
66{
67 const Int_t nrows = row_upb-row_lwb+1;
68 const Int_t ncols = col_upb-col_lwb+1;
69
70 if (nrows < ncols) {
71 Error("TDecompSVD(Int_t,Int_t,Int_t,Int_t","matrix rows should be >= columns");
72 return;
73 }
78 fV.ResizeTo(nrows,ncols); // In the end we only need the nColxnCol part
79}
80
81////////////////////////////////////////////////////////////////////////////////
82/// Constructor for general matrix A .
83
85{
86 R__ASSERT(a.IsValid());
87 if (a.GetNrows() < a.GetNcols()) {
88 Error("TDecompSVD(const TMatrixD &","matrix rows should be >= columns");
89 return;
90 }
91
93 fTol = a.GetTol();
94 if (tol > 0)
95 fTol = tol;
96
97 fRowLwb = a.GetRowLwb();
98 fColLwb = a.GetColLwb();
99 const Int_t nRow = a.GetNrows();
100 const Int_t nCol = a.GetNcols();
101
104 fV.ResizeTo(nRow,nCol); // In the end we only need the nColxnCol part
105
106 fU.UnitMatrix();
107 memcpy(fV.GetMatrixArray(),a.GetMatrixArray(),nRow*nCol*sizeof(Double_t));
108}
109
110////////////////////////////////////////////////////////////////////////////////
111/// Copy constructor
112
117
118////////////////////////////////////////////////////////////////////////////////
119/// SVD decomposition of matrix
120/// If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular
121
123{
124 if (TestBit(kDecomposed)) return kTRUE;
125
126 if ( !TestBit(kMatrixSet) ) {
127 Error("Decompose()","Matrix has not been set");
128 return kFALSE;
129 }
130
131 const Int_t nCol = this->GetNcols();
132 const Int_t rowLwb = this->GetRowLwb();
133 const Int_t colLwb = this->GetColLwb();
134
137 if (nCol > kWorkMax) offDiag.ResizeTo(nCol);
138 else offDiag.Use(nCol,work);
139
140 // step 1: bidiagonalization of A
142 return kFALSE;
143
144 // step 2: diagonalization of bidiagonal matrix
146 return kFALSE;
147
148 // step 3: order singular values and perform permutations
154
155 return kTRUE;
156}
157
158////////////////////////////////////////////////////////////////////////////////
159/// Bidiagonalize the (m x n) - matrix a (stored in v) through a series of Householder
160/// transformations applied to the left (Q^T) and to the right (H) of a ,
161/// so that A = Q . C . H^T with matrix C bidiagonal. Q and H are orthogonal matrices .
162///
163/// Output:
164/// v - (n x n) - matrix H in the (n x n) part of v
165/// u - (m x m) - matrix Q^T
166/// sDiag - diagonal of the (m x n) C
167/// oDiag - off-diagonal elements of matrix C
168///
169/// Test code for the output:
170/// ~~~
171/// const Int_t nRow = v.GetNrows();
172/// const Int_t nCol = v.GetNcols();
173/// TMatrixD H(v); H.ResizeTo(nCol,nCol);
174/// TMatrixD E1(nCol,nCol); E1.UnitMatrix();
175/// TMatrixD Ht(TMatrixDBase::kTransposed,H);
176/// Bool_t ok = kTRUE;
177/// ok &= VerifyMatrixIdentity(Ht * H,E1,kTRUE,1.0e-13);
178/// ok &= VerifyMatrixIdentity(H * Ht,E1,kTRUE,1.0e-13);
179/// TMatrixD E2(nRow,nRow); E2.UnitMatrix();
180/// TMatrixD Qt(u);
181/// TMatrixD Q(TMatrixDBase::kTransposed,Qt);
182/// ok &= VerifyMatrixIdentity(Q * Qt,E2,kTRUE,1.0e-13);
183/// TMatrixD C(nRow,nCol);
184/// TMatrixDDiag(C) = sDiag;
185/// for (Int_t i = 0; i < nCol-1; i++)
186/// C(i,i+1) = oDiag(i+1);
187/// TMatrixD A = Q*C*Ht;
188/// ok &= VerifyMatrixIdentity(A,a,kTRUE,1.0e-13);
189/// ~~~
190
192{
193 const Int_t nRow_v = v.GetNrows();
194 const Int_t nCol_v = v.GetNcols();
195 const Int_t nCol_u = u.GetNcols();
196
199
200 for (Int_t i = 0; i < nCol_v; i++) {
201 // Set up Householder Transformation q(i)
202
203 Double_t up,beta;
206 //if (!DefHouseHolder(vc_i,i,i+1,up,beta))
207 // return kFALSE;
208 DefHouseHolder(vc_i,i,i+1,up,beta);
209
210 // Apply q(i) to v
211 for (Int_t j = i; j < nCol_v; j++) {
213 ApplyHouseHolder(vc_i,up,beta,i,i+1,vc_j);
214 }
215
216 // Apply q(i) to u
217 for (Int_t j = 0; j < nCol_u; j++)
218 {
220 ApplyHouseHolder(vc_i,up,beta,i,i+1,uc_j);
221 }
222 }
223 if (i < nCol_v-2) {
224 // set up Householder Transformation h(i)
225 const TVectorD vr_i = TMatrixDRow_const(v,i);
226
227 //if (!DefHouseHolder(vr_i,i+1,i+2,up,beta))
228 // return kFALSE;
229 DefHouseHolder(vr_i,i+1,i+2,up,beta);
230
231 // save h(i)
232 ups[i] = up;
233 betas[i] = beta;
234
235 // apply h(i) to v
236 for (Int_t j = i; j < nRow_v; j++) {
238 ApplyHouseHolder(vr_i,up,beta,i+1,i+2,vr_j);
239
240 // save elements i+2,...in row j of matrix v
241 if (j == i) {
242 for (Int_t k = i+2; k < nCol_v; k++)
243 vr_j(k) = vr_i(k);
244 }
245 }
246 }
247 }
248
249 // copy diagonal of transformed matrix v to sDiag and upper parallel v to oDiag
250 if (nCol_v > 1) {
251 for (Int_t i = 1; i < nCol_v; i++)
252 oDiag(i) = v(i-1,i);
253 }
254 oDiag(0) = 0.;
256
257 // construct product matrix h = h(1)*h(2)*...*h(nCol_v-1), h(nCol_v-1) = I
258
260 for (Int_t i = nCol_v-1; i >= 0; i--) {
261 if (i < nCol_v-1)
263 TMatrixDRow(v,i).Assign( 0.0 );
264 v(i,i) = 1.;
265
266 if (i < nCol_v-2) {
267 for (Int_t k = i; k < nCol_v; k++) {
268 // householder transformation on k-th column
270 ApplyHouseHolder(vr_i,ups[i],betas[i],i+1,i+2,vc_k);
271 }
272 }
273 }
274
275 return kTRUE;
276}
277
278////////////////////////////////////////////////////////////////////////////////
279/// Diagonalizes in an iterative fashion the bidiagonal matrix C as described through
280/// sDiag and oDiag, so that S' = U'^T . C . V' is diagonal. U' and V' are orthogonal
281/// matrices .
282///
283/// Output:
284/// - v - (n x n) - matrix H . V' in the (n x n) part of v
285/// - u - (m x m) - matrix U'^T . Q^T
286/// - sDiag - diagonal of the (m x n) S'
287///
288/// return convergence flag: 0 -> no convergence
289/// 1 -> convergence
290///
291/// Test code for the output:
292/// ~~~
293/// const Int_t nRow = v.GetNrows();
294/// const Int_t nCol = v.GetNcols();
295/// TMatrixD tmp = v; tmp.ResizeTo(nCol,nCol);
296/// TMatrixD Vprime = Ht*tmp;
297/// TMatrixD Vprimet(TMatrixDBase::kTransposed,Vprime);
298/// TMatrixD Uprimet = u*Q;
299/// TMatrixD Uprime(TMatrixDBase::kTransposed,Uprimet);
300/// TMatrixD Sprime(nRow,nCol);
301/// TMatrixDDiag(Sprime) = sDiag;
302/// ok &= VerifyMatrixIdentity(Uprimet * C * Vprime,Sprime,kTRUE,1.0e-13);
303/// ok &= VerifyMatrixIdentity(Q*Uprime * Sprime * Vprimet * Ht,a,kTRUE,1.0e-13);
304/// ~~~
305
307{
308 Bool_t ok = kTRUE;
309 Int_t niter = 0;
310 Double_t bmx = sDiag(0);
311
312 const Int_t nCol_v = v.GetNcols();
313
314 if (nCol_v > 1) {
315 for (Int_t i = 1; i < nCol_v; i++)
317 }
318
319 const Double_t eps = std::numeric_limits<double>::epsilon();
320
321 const Int_t niterm = 10*nCol_v;
322 for (Int_t k = nCol_v-1; k >= 0; k--) {
323 loop:
324 if (k != 0) {
325 // since sDiag(k) == 0 perform Givens transform with result oDiag[k] = 0
326 if (TMath::Abs(sDiag(k)) < eps*bmx)
327 Diag_1(v,sDiag,oDiag,k);
328
329 // find l (1 <= l <=k) so that either oDiag(l) = 0 or sDiag(l-1) = 0.
330 // In the latter case transform oDiag(l) to zero. In both cases the matrix
331 // splits and the bottom right minor begins with row l.
332 // If no such l is found set l = 1.
333
334 Int_t elzero = 0;
335 Int_t l = 0;
336 for (Int_t ll = k; ll >= 0 ; ll--) {
337 l = ll;
338 if (l == 0) {
339 elzero = 0;
340 break;
341 } else if (TMath::Abs(oDiag(l)) < eps*bmx) {
342 elzero = 1;
343 break;
344 } else if (TMath::Abs(sDiag(l-1)) < eps*bmx)
345 elzero = 0;
346 }
347 if (l > 0 && !elzero)
348 Diag_2(sDiag,oDiag,k,l);
349 if (l != k) {
350 // one more QR pass with order k
351 Diag_3(v,u,sDiag,oDiag,k,l);
352 niter++;
353 if (niter <= niterm) goto loop;
354 ::Error("TDecompSVD::Diagonalize","no convergence after %d steps",niter);
355 ok = kFALSE;
356 }
357 }
358
359 if (sDiag(k) < 0.) {
360 // for negative singular values perform change of sign
361 sDiag(k) = -sDiag(k);
362 TMatrixDColumn(v,k) *= -1.0;
363 }
364 // order is decreased by one in next pass
365 }
366
367 return ok;
368}
369
370////////////////////////////////////////////////////////////////////////////////
371/// Step 1 in the matrix diagonalization
372
374{
375 const Int_t nCol_v = v.GetNcols();
376
378 for (Int_t i = k-1; i >= 0; i--) {
380 Double_t h,cs,sn;
381 if (i == k-1)
382 DefAplGivens(sDiag[i],oDiag[i+1],cs,sn);
383 else
385 if (i > 0) {
386 h = 0.;
387 ApplyGivens(oDiag[i],h,cs,sn);
388 }
389 for (Int_t j = 0; j < nCol_v; j++)
391 }
392}
393
394////////////////////////////////////////////////////////////////////////////////
395/// Step 2 in the matrix diagonalization
396
398{
399 for (Int_t i = l; i <= k; i++) {
400 Double_t h,cs,sn;
401 if (i == l)
403 else
405 if (i < k) {
406 h = 0.;
407 ApplyGivens(oDiag(i+1),h,cs,sn);
408 }
409 }
410}
411
412////////////////////////////////////////////////////////////////////////////////
413/// Step 3 in the matrix diagonalization
414
416{
417 Double_t *pS = sDiag.GetMatrixArray();
418 Double_t *pO = oDiag.GetMatrixArray();
419
420 // determine shift parameter
421
422 const Double_t psk1 = pS[k-1];
423 const Double_t psk = pS[k];
424 const Double_t pok1 = pO[k-1];
425 const Double_t pok = pO[k];
426 const Double_t psl = pS[l];
427
428 Double_t f;
429 if (psl == 0.0 || pok == 0.0 || psk1 == 0.0) {
430 const Double_t b = ((psk1-psk)*(psk1+psk)+pok1*pok1)/2.0;
431 const Double_t c = (psk*pok1)*(psk*pok1);
432
433 Double_t shift = 0.0;
434 if ((b != 0.0) | (c != 0.0)) {
435 shift = TMath::Sqrt(b*b+c);
436 if (b < 0.0)
437 shift = -shift;
438 shift = c/(b+shift);
439 }
440
441 f = (psl+psk)*(psl-psk)+shift;
442 } else {
443 f = ((psk1-psk)*(psk1+psk)+(pok1-pok)*(pok1+pok))/(2.*pok*psk1);
444 const Double_t g = TMath::Hypot(1.,f);
445 const Double_t t = (f >= 0.) ? f+g : f-g;
446
447 f = ((psl-psk)*(psl+psk)+pok*(psk1/t-pok))/psl;
448 }
449
450 const Int_t nCol_v = v.GetNcols();
451 const Int_t nCol_u = u.GetNcols();
452
453 Double_t h,cs,sn;
454 Int_t j;
455 for (Int_t i = l; i <= k-1; i++) {
456 if (i == l)
457 // define r[l]
458 DefGivens(f,pO[i+1],cs,sn);
459 else
460 // define r[i]
461 DefAplGivens(pO[i],h,cs,sn);
462
463 ApplyGivens(pS[i],pO[i+1],cs,sn);
464 h = 0.;
465 ApplyGivens(h,pS[i+1],cs,sn);
466
469 for (j = 0; j < nCol_v; j++)
471 // define t[i]
472 DefAplGivens(pS[i],h,cs,sn);
473 ApplyGivens(pO[i+1],pS[i+1],cs,sn);
474 if (i < k-1) {
475 h = 0.;
476 ApplyGivens(h,pO[i+2],cs,sn);
477 }
478
481 for (j = 0; j < nCol_u; j++)
483 }
484}
485
486////////////////////////////////////////////////////////////////////////////////
487/// Perform a permutation transformation on the diagonal matrix S', so that
488/// matrix S'' = U''^T . S' . V'' has diagonal elements ordered such that they
489/// do not increase.
490///
491/// Output:
492/// - v - (n x n) - matrix H . V' . V'' in the (n x n) part of v
493/// - u - (m x m) - matrix U''^T . U'^T . Q^T
494/// - sDiag - diagonal of the (m x n) S''
495
497{
498 const Int_t nCol_v = v.GetNcols();
499 const Int_t nCol_u = u.GetNcols();
500
501 Double_t *pS = sDiag.GetMatrixArray();
502 Double_t *pV = v.GetMatrixArray();
503 Double_t *pU = u.GetMatrixArray();
504
505 // order singular values
506
507 Int_t i,j;
508 if (nCol_v > 1) {
509 while (true) {
510 Bool_t found = kFALSE;
511 i = 1;
512 while (!found && i < nCol_v) {
513 if (pS[i] > pS[i-1])
514 found = kTRUE;
515 else
516 i++;
517 }
518 if (!found) break;
519 for (i = 1; i < nCol_v; i++) {
520 Double_t t = pS[i-1];
521 Int_t k = i-1;
522 for (j = i; j < nCol_v; j++) {
523 if (t < pS[j]) {
524 t = pS[j];
525 k = j;
526 }
527 }
528 if (k != i-1) {
529 // perform permutation on singular values
530 pS[k] = pS[i-1];
531 pS[i-1] = t;
532 // perform permutation on matrix v
533 for (j = 0; j < nCol_v; j++) {
534 const Int_t off_j = j*nCol_v;
535 t = pV[off_j+k];
536 pV[off_j+k] = pV[off_j+i-1];
537 pV[off_j+i-1] = t;
538 }
539 // perform permutation on vector u
540 for (j = 0; j < nCol_u; j++) {
541 const Int_t off_k = k*nCol_u;
542 const Int_t off_i1 = (i-1)*nCol_u;
543 t = pU[off_k+j];
544 pU[off_k+j] = pU[off_i1+j];
545 pU[off_i1+j] = t;
546 }
547 }
548 }
549 }
550 }
551}
552
553////////////////////////////////////////////////////////////////////////////////
554/// Reconstruct the original matrix using the decomposition parts
555
557{
558 if (TestBit(kSingular)) {
559 Error("GetMatrix()","Matrix is singular");
560 return TMatrixD();
561 }
562 if ( !TestBit(kDecomposed) ) {
563 if (!Decompose()) {
564 Error("GetMatrix()","Decomposition failed");
565 return TMatrixD();
566 }
567 }
568
569 const Int_t nRows = fU.GetNrows();
570 const Int_t nCols = fV.GetNcols();
571 const Int_t colLwb = this->GetColLwb();
575 return fU * s * vt;
576}
577
578////////////////////////////////////////////////////////////////////////////////
579/// Set matrix to be decomposed
580
582{
583 R__ASSERT(a.IsValid());
584
585 ResetStatus();
586 if (a.GetNrows() < a.GetNcols()) {
587 Error("TDecompSVD(const TMatrixD &","matrix rows should be >= columns");
588 return;
589 }
590
592 fCondition = -1.0;
593
594 fRowLwb = a.GetRowLwb();
595 fColLwb = a.GetColLwb();
596 const Int_t nRow = a.GetNrows();
597 const Int_t nCol = a.GetNcols();
598
601 fV.ResizeTo(nRow,nCol); // In the end we only need the nColxnCol part
602
603 fU.UnitMatrix();
604 memcpy(fV.GetMatrixArray(),a.GetMatrixArray(),nRow*nCol*sizeof(Double_t));
605}
606
607////////////////////////////////////////////////////////////////////////////////
608/// Solve Ax=b assuming the SVD form of A is stored . Solution returned in b.
609/// If A is of size (m x n), input vector b should be of size (m), however,
610/// the solution, returned in b, will be in the first (n) elements .
611///
612/// For m > n , x is the least-squares solution of min(A . x - b)
613
615{
616 R__ASSERT(b.IsValid());
617 if (TestBit(kSingular)) {
618 return kFALSE;
619 }
620 if ( !TestBit(kDecomposed) ) {
621 if (!Decompose()) {
622 return kFALSE;
623 }
624 }
625
626 if (fU.GetNrows() != b.GetNrows() || fU.GetRowLwb() != b.GetLwb())
627 {
628 Error("Solve(TVectorD &","vector and matrix incompatible");
629 return kFALSE;
630 }
631
632 // We start with fU fSig fV^T x = b, and turn it into fV^T x = fSig^-1 fU^T b
633 // Form tmp = fSig^-1 fU^T b but ignore diagonal elements with
634 // fSig(i) < fTol * max(fSig)
635
636 const Int_t lwb = fU.GetColLwb();
637 const Int_t upb = lwb+fV.GetNcols()-1;
638 const Double_t threshold = fSig(lwb)*fTol;
639
640 TVectorD tmp(lwb,upb);
641 for (Int_t irow = lwb; irow <= upb; irow++) {
642 Double_t r = 0.0;
643 if (fSig(irow) > threshold) {
645 r = uc_i * b;
646 r /= fSig(irow);
647 }
648 tmp(irow) = r;
649 }
650
651 if (b.GetNrows() > fV.GetNrows()) {
653 tmp2.Use(lwb,upb,b.GetMatrixArray());
654 tmp2 = fV*tmp;
655 } else
656 b = fV*tmp;
657
658 return kTRUE;
659}
660
661////////////////////////////////////////////////////////////////////////////////
662/// Solve Ax=b assuming the SVD form of A is stored . Solution returned in the
663/// matrix column cb b.
664/// If A is of size (m x n), input vector b should be of size (m), however,
665/// the solution, returned in b, will be in the first (n) elements .
666///
667/// For m > n , x is the least-squares solution of min(A . x - b)
668
670{
671 TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
672 R__ASSERT(b->IsValid());
673 if (TestBit(kSingular)) {
674 return kFALSE;
675 }
676 if ( !TestBit(kDecomposed) ) {
677 if (!Decompose()) {
678 return kFALSE;
679 }
680 }
681
682 if (fU.GetNrows() != b->GetNrows() || fU.GetRowLwb() != b->GetRowLwb())
683 {
684 Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
685 return kFALSE;
686 }
687
688 // We start with fU fSig fV^T x = b, and turn it into fV^T x = fSig^-1 fU^T b
689 // Form tmp = fSig^-1 fU^T b but ignore diagonal elements in
690 // fSig(i) < fTol * max(fSig)
691
692 const Int_t lwb = fU.GetColLwb();
693 const Int_t upb = lwb+fV.GetNcols()-1;
694 const Double_t threshold = fSig(lwb)*fTol;
695
696 TVectorD tmp(lwb,upb);
697 const TVectorD vb = cb;
698 for (Int_t irow = lwb; irow <= upb; irow++) {
699 Double_t r = 0.0;
700 if (fSig(irow) > threshold) {
702 r = uc_i * vb;
703 r /= fSig(irow);
704 }
705 tmp(irow) = r;
706 }
707
708 if (b->GetNrows() > fV.GetNrows()) {
709 const TVectorD tmp2 = fV*tmp;
710 TVectorD tmp3(cb);
711 tmp3.SetSub(tmp2.GetLwb(),tmp2);
712 cb = tmp3;
713 } else
714 cb = fV*tmp;
715
716 return kTRUE;
717}
718
719////////////////////////////////////////////////////////////////////////////////
720/// Solve A^T x=b assuming the SVD form of A is stored . Solution returned in b.
721
723{
724 R__ASSERT(b.IsValid());
725 if (TestBit(kSingular)) {
726 return kFALSE;
727 }
728 if ( !TestBit(kDecomposed) ) {
729 if (!Decompose()) {
730 return kFALSE;
731 }
732 }
733
734 if (fU.GetNrows() != fV.GetNrows() || fU.GetRowLwb() != fV.GetRowLwb()) {
735 Error("TransSolve(TVectorD &","matrix should be square");
736 return kFALSE;
737 }
738
739 if (fV.GetNrows() != b.GetNrows() || fV.GetRowLwb() != b.GetLwb())
740 {
741 Error("TransSolve(TVectorD &","vector and matrix incompatible");
742 return kFALSE;
743 }
744
745 // We start with fV fSig fU^T x = b, and turn it into fU^T x = fSig^-1 fV^T b
746 // Form tmp = fSig^-1 fV^T b but ignore diagonal elements in
747 // fSig(i) < fTol * max(fSig)
748
749 const Int_t lwb = fU.GetColLwb();
750 const Int_t upb = lwb+fV.GetNcols()-1;
751 const Double_t threshold = fSig(lwb)*fTol;
752
753 TVectorD tmp(lwb,upb);
754 for (Int_t i = lwb; i <= upb; i++) {
755 Double_t r = 0.0;
756 if (fSig(i) > threshold) {
757 const TVectorD vc = TMatrixDColumn(fV,i);
758 r = vc * b;
759 r /= fSig(i);
760 }
761 tmp(i) = r;
762 }
763 b = fU*tmp;
764
765 return kTRUE;
766}
767
768////////////////////////////////////////////////////////////////////////////////
769/// Solve A^T x=b assuming the SVD form of A is stored . Solution returned in b.
770
772{
773 TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
774 R__ASSERT(b->IsValid());
775 if (TestBit(kSingular)) {
776 return kFALSE;
777 }
778 if ( !TestBit(kDecomposed) ) {
779 if (!Decompose()) {
780 return kFALSE;
781 }
782 }
783
784 if (fU.GetNrows() != fV.GetNrows() || fU.GetRowLwb() != fV.GetRowLwb()) {
785 Error("TransSolve(TMatrixDColumn &","matrix should be square");
786 return kFALSE;
787 }
788
789 if (fV.GetNrows() != b->GetNrows() || fV.GetRowLwb() != b->GetRowLwb())
790 {
791 Error("TransSolve(TMatrixDColumn &","vector and matrix incompatible");
792 return kFALSE;
793 }
794
795 // We start with fV fSig fU^T x = b, and turn it into fU^T x = fSig^-1 fV^T b
796 // Form tmp = fSig^-1 fV^T b but ignore diagonal elements in
797 // fSig(i) < fTol * max(fSig)
798
799 const Int_t lwb = fU.GetColLwb();
800 const Int_t upb = lwb+fV.GetNcols()-1;
801 const Double_t threshold = fSig(lwb)*fTol;
802
803 const TVectorD vb = cb;
804 TVectorD tmp(lwb,upb);
805 for (Int_t i = lwb; i <= upb; i++) {
806 Double_t r = 0.0;
807 if (fSig(i) > threshold) {
808 const TVectorD vc = TMatrixDColumn(fV,i);
809 r = vc * vb;
810 r /= fSig(i);
811 }
812 tmp(i) = r;
813 }
814 cb = fU*tmp;
815
816 return kTRUE;
817}
818
819////////////////////////////////////////////////////////////////////////////////
820/// Matrix condition number
821
823{
824 if ( !TestBit(kCondition) ) {
825 fCondition = -1;
826 if (TestBit(kSingular))
827 return fCondition;
828 if ( !TestBit(kDecomposed) ) {
829 if (!Decompose())
830 return fCondition;
831 }
832 const Int_t colLwb = GetColLwb();
833 const Int_t nCols = GetNcols();
834 const Double_t max = fSig(colLwb);
835 const Double_t min = fSig(colLwb+nCols-1);
836 fCondition = (min > 0.0) ? max/min : -1.0;
838 }
839 return fCondition;
840}
841
842////////////////////////////////////////////////////////////////////////////////
843/// Matrix determinant det = d1*TMath::Power(2.,d2)
844
846{
847 if ( !TestBit(kDetermined) ) {
848 if ( !TestBit(kDecomposed) )
849 Decompose();
850 if (TestBit(kSingular)) {
851 fDet1 = 0.0;
852 fDet2 = 0.0;
853 } else {
855 }
857 }
858 d1 = fDet1;
859 d2 = fDet2;
860}
861
862////////////////////////////////////////////////////////////////////////////////
863
865{
866 return fU.GetNrows();
867}
868
870{
871 return fV.GetNcols();
872}
873
874////////////////////////////////////////////////////////////////////////////////
875/// For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
876/// The user should always supply a matrix of size (m x m) !
877/// If m > n , only the (n x m) part of the returned (pseudo inverse) matrix
878/// should be used .
879
881{
882 const Int_t rowLwb = GetRowLwb();
883 const Int_t colLwb = GetColLwb();
884 const Int_t nRows = fU.GetNrows();
885
886 if (inv.GetNrows() != nRows || inv.GetNcols() != nRows ||
887 inv.GetRowLwb() != rowLwb || inv.GetColLwb() != colLwb) {
888 Error("Invert(TMatrixD &","Input matrix has wrong shape");
889 return kFALSE;
890 }
891
892 inv.UnitMatrix();
893 Bool_t status = MultiSolve(inv);
894
895 return status;
896}
897
898////////////////////////////////////////////////////////////////////////////////
899/// For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
900/// (n x m) Ainv is returned .
901
903{
904 const Int_t rowLwb = GetRowLwb();
905 const Int_t colLwb = GetColLwb();
906 const Int_t rowUpb = rowLwb+fU.GetNrows()-1;
908 inv.UnitMatrix();
909 status = MultiSolve(inv);
910 inv.ResizeTo(rowLwb,rowLwb+fV.GetNcols()-1,colLwb,colLwb+fU.GetNrows()-1);
911
912 return inv;
913}
914
915////////////////////////////////////////////////////////////////////////////////
916/// Print class members
917
919{
921 fU.Print("fU");
922 fV.Print("fV");
923 fSig.Print("fSig");
924}
925
926////////////////////////////////////////////////////////////////////////////////
927/// Assignment operator
928
930{
931 if (this != &source) {
933 fU.ResizeTo(source.fU);
934 fU = source.fU;
935 fV.ResizeTo(source.fV);
936 fV = source.fV;
937 fSig.ResizeTo(source.fSig);
938 fSig = source.fSig;
939 }
940 return *this;
941}
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Bool_t DefHouseHolder(const TVectorD &vc, Int_t lp, Int_t l, Double_t &up, Double_t &b, Double_t tol=0.0)
Define a Householder-transformation through the parameters up and b .
void ApplyGivens(Double_t &z1, Double_t &z2, Double_t c, Double_t s)
Apply a Givens transformation as defined by c and s to the vector components v1 and v2 .
void ApplyHouseHolder(const TVectorD &vc, Double_t up, Double_t b, Int_t lp, Int_t l, TMatrixDRow &cr)
Apply Householder-transformation.
void DefGivens(Double_t v1, Double_t v2, Double_t &c, Double_t &s)
Defines a Givens-rotation by calculating 2 rotation parameters c and s.
void DefAplGivens(Double_t &v1, Double_t &v2, Double_t &c, Double_t &s)
Define and apply a Givens-rotation by calculating 2 rotation parameters c and s.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
TMatrixTRow< Double_t > TMatrixDRow
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
TMatrixTDiag< Double_t > TMatrixDDiag
TMatrixTRow_const< Double_t > TMatrixDRow_const
TMatrixTColumn< Double_t > TMatrixDColumn
TMatrixT< Double_t > TMatrixD
Definition TMatrixDfwd.h:23
Array of doubles (64 bits per element).
Definition TArrayD.h:27
Decomposition Base class.
Definition TDecompBase.h:34
static void DiagProd(const TVectorD &diag, Double_t tol, Double_t &d1, Double_t &d2)
Int_t GetRowLwb() const
Definition TDecompBase.h:73
Double_t fDet1
Definition TDecompBase.h:37
Double_t fDet2
Definition TDecompBase.h:38
void ResetStatus()
Definition TDecompBase.h:43
Int_t GetColLwb() const
Definition TDecompBase.h:74
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
Int_t fRowLwb
Definition TDecompBase.h:40
Double_t fTol
Definition TDecompBase.h:36
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Double_t fCondition
Definition TDecompBase.h:39
Int_t fColLwb
Definition TDecompBase.h:41
void Print(Option_t *opt="") const override
Print class members.
Single Value Decomposition class.
Definition TDecompSVD.h:24
void Det(Double_t &d1, Double_t &d2) override
Matrix determinant det = d1*TMath::Power(2.,d2)
TVectorD fSig
Definition TDecompSVD.h:30
static void Diag_3(TMatrixD &v, TMatrixD &u, TVectorD &sDiag, TVectorD &oDiag, Int_t k, Int_t l)
Step 3 in the matrix diagonalization.
static void Diag_2(TVectorD &sDiag, TVectorD &oDiag, Int_t k, Int_t l)
Step 2 in the matrix diagonalization.
TDecompSVD & operator=(const TDecompSVD &source)
Assignment operator.
Bool_t Solve(TVectorD &b) override
Solve Ax=b assuming the SVD form of A is stored .
TMatrixD Invert()
Definition TDecompSVD.h:82
TMatrixD fV
Definition TDecompSVD.h:29
static Bool_t Bidiagonalize(TMatrixD &v, TMatrixD &u, TVectorD &sDiag, TVectorD &oDiag)
Bidiagonalize the (m x n) - matrix a (stored in v) through a series of Householder transformations ap...
static void Diag_1(TMatrixD &v, TVectorD &sDiag, TVectorD &oDiag, Int_t k)
Step 1 in the matrix diagonalization.
TMatrixD fU
Definition TDecompSVD.h:28
static Bool_t Diagonalize(TMatrixD &v, TMatrixD &u, TVectorD &sDiag, TVectorD &oDiag)
Diagonalizes in an iterative fashion the bidiagonal matrix C as described through sDiag and oDiag,...
Double_t Condition() override
Matrix condition number.
Bool_t Decompose() override
SVD decomposition of matrix If the decomposition succeeds, bit kDecomposed is set ,...
Int_t GetNrows() const override
const TMatrixD GetMatrix()
Reconstruct the original matrix using the decomposition parts.
Bool_t TransSolve(TVectorD &b) override
Solve A^T x=b assuming the SVD form of A is stored . Solution returned in b.
virtual void SetMatrix(const TMatrixD &a)
Set matrix to be decomposed.
Int_t GetNcols() const override
static void SortSingular(TMatrixD &v, TMatrixD &u, TVectorD &sDiag)
Perform a permutation transformation on the diagonal matrix S', so that matrix S'' = U''^T .
void Print(Option_t *opt="") const override
Print class members.
Int_t GetNrows() const
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Int_t GetRowLwb() const
Int_t GetColLwb() const
Int_t GetNcols() const
virtual TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift)
Shift the row index by adding row_shift and the column index by adding col_shift, respectively.
const TMatrixTBase< Element > * GetMatrix() const
const Element * GetMatrixArray() const override
Definition TMatrixT.h:228
TMatrixT< Element > & Transpose(const TMatrixT< Element > &source)
Transpose matrix source.
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:202
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:864
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
TVectorT< Element > & Shift(Int_t row_shift)
Definition TVectorT.h:107
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:293
void Print(Option_t *option="") const override
Print the vector as a list of elements.
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Double_t Hypot(Double_t x, Double_t y)
Returns sqrt(x*x + y*y)
Definition TMath.cxx:59
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition rsaaux.cxx:949
TLine l
Definition textangle.C:4