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