Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 ClassImp(TDecompSVD);
48 
49 ////////////////////////////////////////////////////////////////////////////////
50 /// Constructor for (nrows x ncols) matrix
51 
52 TDecompSVD::TDecompSVD(Int_t nrows,Int_t ncols)
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 
66 TDecompSVD::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 
85 TDecompSVD::TDecompSVD(const TMatrixD &a,Double_t tol)
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 
93  SetBit(kMatrixSet);
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 
114 TDecompSVD::TDecompSVD(const TDecompSVD &another): TDecompBase(another)
115 {
116  *this = another;
117 }
118 
119 ////////////////////////////////////////////////////////////////////////////////
120 /// SVD decomposition of matrix
121 /// If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular
122 
123 Bool_t TDecompSVD::Decompose()
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
150  SortSingular(fV,fU,fSig);
151  fV.ResizeTo(nCol,nCol); fV.Shift(colLwb,colLwb);
152  fSig.Shift(colLwb);
153  fU.Transpose(fU); fU.Shift(rowLwb,colLwb);
154  SetBit(kDecomposed);
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 
192 Bool_t TDecompSVD::Bidiagonalize(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag)
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++) {
213  TMatrixDColumn vc_j = TMatrixDColumn(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
270  TMatrixDColumn vc_k = TMatrixDColumn(v,k);
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 
307 Bool_t TDecompSVD::Diagonalize(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag)
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 
374 void TDecompSVD::Diag_1(TMatrixD &v,TVectorD &sDiag,TVectorD &oDiag,Int_t k)
375 {
376  const Int_t nCol_v = v.GetNcols();
377 
378  TMatrixDColumn vc_k = TMatrixDColumn(v,k);
379  for (Int_t i = k-1; i >= 0; i--) {
380  TMatrixDColumn vc_i = TMatrixDColumn(v,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 
398 void TDecompSVD::Diag_2(TVectorD &sDiag,TVectorD &oDiag,Int_t k,Int_t l)
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 
416 void TDecompSVD::Diag_3(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag,Int_t k,Int_t l)
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 
468  TMatrixDColumn vc_i = TMatrixDColumn(v,i);
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 
497 void TDecompSVD::SortSingular(TMatrixD &v,TMatrixD &u,TVectorD &sDiag)
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 
557 const TMatrixD TDecompSVD::GetMatrix()
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;
575  const TMatrixD vt(TMatrixD::kTransposed,fV);
576  return fU * s * vt;
577 }
578 
579 ////////////////////////////////////////////////////////////////////////////////
580 /// Set matrix to be decomposed
581 
582 void TDecompSVD::SetMatrix(const TMatrixD &a)
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 
592  SetBit(kMatrixSet);
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 
615 Bool_t TDecompSVD::Solve(TVectorD &b)
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 
670 Bool_t TDecompSVD::Solve(TMatrixDColumn &cb)
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 
723 Bool_t TDecompSVD::TransSolve(TVectorD &b)
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 
772 Bool_t TDecompSVD::TransSolve(TMatrixDColumn &cb)
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 
823 Double_t TDecompSVD::Condition()
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;
838  SetBit(kCondition);
839  }
840  return fCondition;
841 }
842 
843 ////////////////////////////////////////////////////////////////////////////////
844 /// Matrix determinant det = d1*TMath::Power(2.,d2)
845 
846 void TDecompSVD::Det(Double_t &d1,Double_t &d2)
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 {
855  DiagProd(fSig,fTol,fDet1,fDet2);
856  }
857  SetBit(kDetermined);
858  }
859  d1 = fDet1;
860  d2 = fDet2;
861 }
862 
863 ////////////////////////////////////////////////////////////////////////////////
864 
865 Int_t TDecompSVD::GetNrows () const
866 {
867  return fU.GetNrows();
868 }
869 
870 Int_t TDecompSVD::GetNcols () const
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 
881 Bool_t TDecompSVD::Invert(TMatrixD &inv)
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 
903 TMatrixD TDecompSVD::Invert(Bool_t &status)
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 
919 void TDecompSVD::Print(Option_t *opt) const
920 {
921  TDecompBase::Print(opt);
922  fU.Print("fU");
923  fV.Print("fV");
924  fSig.Print("fSig");
925 }
926 
927 ////////////////////////////////////////////////////////////////////////////////
928 /// Assignment operator
929 
930 TDecompSVD &TDecompSVD::operator=(const TDecompSVD &source)
931 {
932  if (this != &source) {
933  TDecompBase::operator=(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 }