Eigen  3.3.90 (mercurial changeset 493691b29be1)
BDCSVD.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // We used the "A Divide-And-Conquer Algorithm for the Bidiagonal SVD"
5 // research report written by Ming Gu and Stanley C.Eisenstat
6 // The code variable names correspond to the names they used in their
7 // report
8 //
9 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
10 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
11 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
12 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
13 // Copyright (C) 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
14 // Copyright (C) 2014-2017 Gael Guennebaud <gael.guennebaud@inria.fr>
15 //
16 // Source Code Form is subject to the terms of the Mozilla
17 // Public License v. 2.0. If a copy of the MPL was not distributed
18 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
19 
20 #ifndef EIGEN_BDCSVD_H
21 #define EIGEN_BDCSVD_H
22 // #define EIGEN_BDCSVD_DEBUG_VERBOSE
23 // #define EIGEN_BDCSVD_SANITY_CHECKS
24 
25 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
26 #undef eigen_internal_assert
27 #define eigen_internal_assert(X) assert(X);
28 #endif
29 
30 namespace Eigen {
31 
32 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
33 IOFormat bdcsvdfmt(8, 0, ", ", "\n", " [", "]");
34 #endif
35 
36 template<typename _MatrixType> class BDCSVD;
37 
38 namespace internal {
39 
40 template<typename _MatrixType>
41 struct traits<BDCSVD<_MatrixType> >
42 {
43  typedef _MatrixType MatrixType;
44 };
45 
46 } // end namespace internal
47 
48 
71 template<typename _MatrixType>
72 class BDCSVD : public SVDBase<BDCSVD<_MatrixType> >
73 {
74  typedef SVDBase<BDCSVD> Base;
75 
76 public:
77  using Base::rows;
78  using Base::cols;
79  using Base::computeU;
80  using Base::computeV;
81 
82  typedef _MatrixType MatrixType;
83  typedef typename MatrixType::Scalar Scalar;
84  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
85  typedef typename NumTraits<RealScalar>::Literal Literal;
86  enum {
87  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
88  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
89  DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime),
90  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
91  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
92  MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime, MaxColsAtCompileTime),
93  MatrixOptions = MatrixType::Options
94  };
95 
96  typedef typename Base::MatrixUType MatrixUType;
97  typedef typename Base::MatrixVType MatrixVType;
98  typedef typename Base::SingularValuesType SingularValuesType;
99 
100  typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> MatrixX;
101  typedef Matrix<RealScalar, Dynamic, Dynamic, ColMajor> MatrixXr;
102  typedef Matrix<RealScalar, Dynamic, 1> VectorType;
103  typedef Array<RealScalar, Dynamic, 1> ArrayXr;
104  typedef Array<Index,1,Dynamic> ArrayXi;
105  typedef Ref<ArrayXr> ArrayRef;
106  typedef Ref<ArrayXi> IndicesRef;
107 
113  BDCSVD() : m_algoswap(16), m_numIters(0)
114  {}
115 
116 
123  BDCSVD(Index rows, Index cols, unsigned int computationOptions = 0)
124  : m_algoswap(16), m_numIters(0)
125  {
126  allocate(rows, cols, computationOptions);
127  }
128 
139  BDCSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
140  : m_algoswap(16), m_numIters(0)
141  {
142  compute(matrix, computationOptions);
143  }
144 
145  ~BDCSVD()
146  {
147  }
148 
159  BDCSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
160 
167  BDCSVD& compute(const MatrixType& matrix)
168  {
169  return compute(matrix, this->m_computationOptions);
170  }
171 
172  void setSwitchSize(int s)
173  {
174  eigen_assert(s>3 && "BDCSVD the size of the algo switch has to be greater than 3");
175  m_algoswap = s;
176  }
177 
178 private:
179  void allocate(Index rows, Index cols, unsigned int computationOptions);
180  void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift);
181  void computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V);
182  void computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, VectorType& singVals, ArrayRef shifts, ArrayRef mus);
183  void perturbCol0(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, const VectorType& singVals, const ArrayRef& shifts, const ArrayRef& mus, ArrayRef zhat);
184  void computeSingVecs(const ArrayRef& zhat, const ArrayRef& diag, const IndicesRef& perm, const VectorType& singVals, const ArrayRef& shifts, const ArrayRef& mus, MatrixXr& U, MatrixXr& V);
185  void deflation43(Index firstCol, Index shift, Index i, Index size);
186  void deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size);
187  void deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift);
188  template<typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV>
189  void copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naivev);
190  void structured_update(Block<MatrixXr,Dynamic,Dynamic> A, const MatrixXr &B, Index n1);
191  static RealScalar secularEq(RealScalar x, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const ArrayRef& diagShifted, RealScalar shift);
192 
193 protected:
194  MatrixXr m_naiveU, m_naiveV;
195  MatrixXr m_computed;
196  Index m_nRec;
197  ArrayXr m_workspace;
198  ArrayXi m_workspaceI;
199  int m_algoswap;
200  bool m_isTranspose, m_compU, m_compV;
201 
202  using Base::m_singularValues;
203  using Base::m_diagSize;
204  using Base::m_computeFullU;
205  using Base::m_computeFullV;
206  using Base::m_computeThinU;
207  using Base::m_computeThinV;
208  using Base::m_matrixU;
209  using Base::m_matrixV;
210  using Base::m_isInitialized;
211  using Base::m_nonzeroSingularValues;
212 
213 public:
214  int m_numIters;
215 }; //end class BDCSVD
216 
217 
218 // Method to allocate and initialize matrix and attributes
219 template<typename MatrixType>
220 void BDCSVD<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
221 {
222  m_isTranspose = (cols > rows);
223 
224  if (Base::allocate(rows, cols, computationOptions))
225  return;
226 
227  m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize );
228  m_compU = computeV();
229  m_compV = computeU();
230  if (m_isTranspose)
231  std::swap(m_compU, m_compV);
232 
233  if (m_compU) m_naiveU = MatrixXr::Zero(m_diagSize + 1, m_diagSize + 1 );
234  else m_naiveU = MatrixXr::Zero(2, m_diagSize + 1 );
235 
236  if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize);
237 
238  m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3);
239  m_workspaceI.resize(3*m_diagSize);
240 }// end allocate
241 
242 template<typename MatrixType>
243 BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions)
244 {
245 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
246  std::cout << "\n\n\n======================================================================================================================\n\n\n";
247 #endif
248  allocate(matrix.rows(), matrix.cols(), computationOptions);
249  using std::abs;
250 
251  const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
252 
253  //**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return
254  if(matrix.cols() < m_algoswap)
255  {
256  // FIXME this line involves temporaries
257  JacobiSVD<MatrixType> jsvd(matrix,computationOptions);
258  if(computeU()) m_matrixU = jsvd.matrixU();
259  if(computeV()) m_matrixV = jsvd.matrixV();
260  m_singularValues = jsvd.singularValues();
261  m_nonzeroSingularValues = jsvd.nonzeroSingularValues();
262  m_isInitialized = true;
263  return *this;
264  }
265 
266  //**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows
267  RealScalar scale = matrix.cwiseAbs().maxCoeff();
268  if(scale==Literal(0)) scale = Literal(1);
269  MatrixX copy;
270  if (m_isTranspose) copy = matrix.adjoint()/scale;
271  else copy = matrix/scale;
272 
273  //**** step 1 - Bidiagonalization
274  // FIXME this line involves temporaries
275  internal::UpperBidiagonalization<MatrixX> bid(copy);
276 
277  //**** step 2 - Divide & Conquer
278  m_naiveU.setZero();
279  m_naiveV.setZero();
280  // FIXME this line involves a temporary matrix
281  m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
282  m_computed.template bottomRows<1>().setZero();
283  divide(0, m_diagSize - 1, 0, 0, 0);
284 
285  //**** step 3 - Copy singular values and vectors
286  for (int i=0; i<m_diagSize; i++)
287  {
288  RealScalar a = abs(m_computed.coeff(i, i));
289  m_singularValues.coeffRef(i) = a * scale;
290  if (a<considerZero)
291  {
292  m_nonzeroSingularValues = i;
293  m_singularValues.tail(m_diagSize - i - 1).setZero();
294  break;
295  }
296  else if (i == m_diagSize - 1)
297  {
298  m_nonzeroSingularValues = i + 1;
299  break;
300  }
301  }
302 
303 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
304 // std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
305 // std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
306 #endif
307  if(m_isTranspose) copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU);
308  else copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV);
309 
310  m_isInitialized = true;
311  return *this;
312 }// end compute
313 
314 
315 template<typename MatrixType>
316 template<typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV>
317 void BDCSVD<MatrixType>::copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naiveV)
318 {
319  // Note exchange of U and V: m_matrixU is set from m_naiveV and vice versa
320  if (computeU())
321  {
322  Index Ucols = m_computeThinU ? m_diagSize : householderU.cols();
323  m_matrixU = MatrixX::Identity(householderU.cols(), Ucols);
324  m_matrixU.topLeftCorner(m_diagSize, m_diagSize) = naiveV.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
325  householderU.applyThisOnTheLeft(m_matrixU); // FIXME this line involves a temporary buffer
326  }
327  if (computeV())
328  {
329  Index Vcols = m_computeThinV ? m_diagSize : householderV.cols();
330  m_matrixV = MatrixX::Identity(householderV.cols(), Vcols);
331  m_matrixV.topLeftCorner(m_diagSize, m_diagSize) = naiveU.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
332  householderV.applyThisOnTheLeft(m_matrixV); // FIXME this line involves a temporary buffer
333  }
334 }
335 
344 template<typename MatrixType>
346 {
347  Index n = A.rows();
348  if(n>100)
349  {
350  // If the matrices are large enough, let's exploit the sparse structure of A by
351  // splitting it in half (wrt n1), and packing the non-zero columns.
352  Index n2 = n - n1;
353  Map<MatrixXr> A1(m_workspace.data() , n1, n);
354  Map<MatrixXr> A2(m_workspace.data()+ n1*n, n2, n);
355  Map<MatrixXr> B1(m_workspace.data()+ n*n, n, n);
356  Map<MatrixXr> B2(m_workspace.data()+2*n*n, n, n);
357  Index k1=0, k2=0;
358  for(Index j=0; j<n; ++j)
359  {
360  if( (A.col(j).head(n1).array()!=Literal(0)).any() )
361  {
362  A1.col(k1) = A.col(j).head(n1);
363  B1.row(k1) = B.row(j);
364  ++k1;
365  }
366  if( (A.col(j).tail(n2).array()!=Literal(0)).any() )
367  {
368  A2.col(k2) = A.col(j).tail(n2);
369  B2.row(k2) = B.row(j);
370  ++k2;
371  }
372  }
373 
374  A.topRows(n1).noalias() = A1.leftCols(k1) * B1.topRows(k1);
375  A.bottomRows(n2).noalias() = A2.leftCols(k2) * B2.topRows(k2);
376  }
377  else
378  {
379  Map<MatrixXr,Aligned> tmp(m_workspace.data(),n,n);
380  tmp.noalias() = A*B;
381  A = tmp;
382  }
383 }
384 
385 // The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the
386 // place of the submatrix we are currently working on.
387 
388 //@param firstCol : The Index of the first column of the submatrix of m_computed and for m_naiveU;
389 //@param lastCol : The Index of the last column of the submatrix of m_computed and for m_naiveU;
390 // lastCol + 1 - firstCol is the size of the submatrix.
391 //@param firstRowW : The Index of the first row of the matrix W that we are to change. (see the reference paper section 1 for more information on W)
392 //@param firstRowW : Same as firstRowW with the column.
393 //@param shift : Each time one takes the left submatrix, one must add 1 to the shift. Why? Because! We actually want the last column of the U submatrix
394 // to become the first column (*coeff) and to shift all the other columns to the right. There are more details on the reference paper.
395 template<typename MatrixType>
396 void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift)
397 {
398  // requires rows = cols + 1;
399  using std::pow;
400  using std::sqrt;
401  using std::abs;
402  const Index n = lastCol - firstCol + 1;
403  const Index k = n/2;
404  const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
405  RealScalar alphaK;
406  RealScalar betaK;
407  RealScalar r0;
408  RealScalar lambda, phi, c0, s0;
409  VectorType l, f;
410  // We use the other algorithm which is more efficient for small
411  // matrices.
412  if (n < m_algoswap)
413  {
414  // FIXME this line involves temporaries
415  JacobiSVD<MatrixXr> b(m_computed.block(firstCol, firstCol, n + 1, n), ComputeFullU | (m_compV ? ComputeFullV : 0));
416  if (m_compU)
417  m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.matrixU();
418  else
419  {
420  m_naiveU.row(0).segment(firstCol, n + 1).real() = b.matrixU().row(0);
421  m_naiveU.row(1).segment(firstCol, n + 1).real() = b.matrixU().row(n);
422  }
423  if (m_compV) m_naiveV.block(firstRowW, firstColW, n, n).real() = b.matrixV();
424  m_computed.block(firstCol + shift, firstCol + shift, n + 1, n).setZero();
425  m_computed.diagonal().segment(firstCol + shift, n) = b.singularValues().head(n);
426  return;
427  }
428  // We use the divide and conquer algorithm
429  alphaK = m_computed(firstCol + k, firstCol + k);
430  betaK = m_computed(firstCol + k + 1, firstCol + k);
431  // The divide must be done in that order in order to have good results. Divide change the data inside the submatrices
432  // and the divide of the right submatrice reads one column of the left submatrice. That's why we need to treat the
433  // right submatrix before the left one.
434  divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift);
435  divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
436 
437  if (m_compU)
438  {
439  lambda = m_naiveU(firstCol + k, firstCol + k);
440  phi = m_naiveU(firstCol + k + 1, lastCol + 1);
441  }
442  else
443  {
444  lambda = m_naiveU(1, firstCol + k);
445  phi = m_naiveU(0, lastCol + 1);
446  }
447  r0 = sqrt((abs(alphaK * lambda) * abs(alphaK * lambda)) + abs(betaK * phi) * abs(betaK * phi));
448  if (m_compU)
449  {
450  l = m_naiveU.row(firstCol + k).segment(firstCol, k);
451  f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1);
452  }
453  else
454  {
455  l = m_naiveU.row(1).segment(firstCol, k);
456  f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1);
457  }
458  if (m_compV) m_naiveV(firstRowW+k, firstColW) = Literal(1);
459  if (r0<considerZero)
460  {
461  c0 = Literal(1);
462  s0 = Literal(0);
463  }
464  else
465  {
466  c0 = alphaK * lambda / r0;
467  s0 = betaK * phi / r0;
468  }
469 
470 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
471  assert(m_naiveU.allFinite());
472  assert(m_naiveV.allFinite());
473  assert(m_computed.allFinite());
474 #endif
475 
476  if (m_compU)
477  {
478  MatrixXr q1 (m_naiveU.col(firstCol + k).segment(firstCol, k + 1));
479  // we shiftW Q1 to the right
480  for (Index i = firstCol + k - 1; i >= firstCol; i--)
481  m_naiveU.col(i + 1).segment(firstCol, k + 1) = m_naiveU.col(i).segment(firstCol, k + 1);
482  // we shift q1 at the left with a factor c0
483  m_naiveU.col(firstCol).segment( firstCol, k + 1) = (q1 * c0);
484  // last column = q1 * - s0
485  m_naiveU.col(lastCol + 1).segment(firstCol, k + 1) = (q1 * ( - s0));
486  // first column = q2 * s0
487  m_naiveU.col(firstCol).segment(firstCol + k + 1, n - k) = m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) * s0;
488  // q2 *= c0
489  m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) *= c0;
490  }
491  else
492  {
493  RealScalar q1 = m_naiveU(0, firstCol + k);
494  // we shift Q1 to the right
495  for (Index i = firstCol + k - 1; i >= firstCol; i--)
496  m_naiveU(0, i + 1) = m_naiveU(0, i);
497  // we shift q1 at the left with a factor c0
498  m_naiveU(0, firstCol) = (q1 * c0);
499  // last column = q1 * - s0
500  m_naiveU(0, lastCol + 1) = (q1 * ( - s0));
501  // first column = q2 * s0
502  m_naiveU(1, firstCol) = m_naiveU(1, lastCol + 1) *s0;
503  // q2 *= c0
504  m_naiveU(1, lastCol + 1) *= c0;
505  m_naiveU.row(1).segment(firstCol + 1, k).setZero();
506  m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).setZero();
507  }
508 
509 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
510  assert(m_naiveU.allFinite());
511  assert(m_naiveV.allFinite());
512  assert(m_computed.allFinite());
513 #endif
514 
515  m_computed(firstCol + shift, firstCol + shift) = r0;
516  m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) = alphaK * l.transpose().real();
517  m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) = betaK * f.transpose().real();
518 
519 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
520  ArrayXr tmp1 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
521 #endif
522  // Second part: try to deflate singular values in combined matrix
523  deflation(firstCol, lastCol, k, firstRowW, firstColW, shift);
524 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
525  ArrayXr tmp2 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
526  std::cout << "\n\nj1 = " << tmp1.transpose().format(bdcsvdfmt) << "\n";
527  std::cout << "j2 = " << tmp2.transpose().format(bdcsvdfmt) << "\n\n";
528  std::cout << "err: " << ((tmp1-tmp2).abs()>1e-12*tmp2.abs()).transpose() << "\n";
529  static int count = 0;
530  std::cout << "# " << ++count << "\n\n";
531  assert((tmp1-tmp2).matrix().norm() < 1e-14*tmp2.matrix().norm());
532 // assert(count<681);
533 // assert(((tmp1-tmp2).abs()<1e-13*tmp2.abs()).all());
534 #endif
535 
536  // Third part: compute SVD of combined matrix
537  MatrixXr UofSVD, VofSVD;
538  VectorType singVals;
539  computeSVDofM(firstCol + shift, n, UofSVD, singVals, VofSVD);
540 
541 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
542  assert(UofSVD.allFinite());
543  assert(VofSVD.allFinite());
544 #endif
545 
546  if (m_compU)
547  structured_update(m_naiveU.block(firstCol, firstCol, n + 1, n + 1), UofSVD, (n+2)/2);
548  else
549  {
550  Map<Matrix<RealScalar,2,Dynamic>,Aligned> tmp(m_workspace.data(),2,n+1);
551  tmp.noalias() = m_naiveU.middleCols(firstCol, n+1) * UofSVD;
552  m_naiveU.middleCols(firstCol, n + 1) = tmp;
553  }
554 
555  if (m_compV) structured_update(m_naiveV.block(firstRowW, firstColW, n, n), VofSVD, (n+1)/2);
556 
557 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
558  assert(m_naiveU.allFinite());
559  assert(m_naiveV.allFinite());
560  assert(m_computed.allFinite());
561 #endif
562 
563  m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero();
564  m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals;
565 }// end divide
566 
567 // Compute SVD of m_computed.block(firstCol, firstCol, n + 1, n); this block only has non-zeros in
568 // the first column and on the diagonal and has undergone deflation, so diagonal is in increasing
569 // order except for possibly the (0,0) entry. The computed SVD is stored U, singVals and V, except
570 // that if m_compV is false, then V is not computed. Singular values are sorted in decreasing order.
571 //
572 // TODO Opportunities for optimization: better root finding algo, better stopping criterion, better
573 // handling of round-off errors, be consistent in ordering
574 // For instance, to solve the secular equation using FMM, see http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
575 template <typename MatrixType>
576 void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V)
577 {
578  const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
579  using std::abs;
580  ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
581  m_workspace.head(n) = m_computed.block(firstCol, firstCol, n, n).diagonal();
582  ArrayRef diag = m_workspace.head(n);
583  diag(0) = Literal(0);
584 
585  // Allocate space for singular values and vectors
586  singVals.resize(n);
587  U.resize(n+1, n+1);
588  if (m_compV) V.resize(n, n);
589 
590 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
591  if (col0.hasNaN() || diag.hasNaN())
592  std::cout << "\n\nHAS NAN\n\n";
593 #endif
594 
595  // Many singular values might have been deflated, the zero ones have been moved to the end,
596  // but others are interleaved and we must ignore them at this stage.
597  // To this end, let's compute a permutation skipping them:
598  Index actual_n = n;
599  while(actual_n>1 && diag(actual_n-1)==Literal(0)) {--actual_n; eigen_internal_assert(col0(actual_n)==Literal(0)); }
600  Index m = 0; // size of the deflated problem
601  for(Index k=0;k<actual_n;++k)
602  if(abs(col0(k))>considerZero)
603  m_workspaceI(m++) = k;
604  Map<ArrayXi> perm(m_workspaceI.data(),m);
605 
606  Map<ArrayXr> shifts(m_workspace.data()+1*n, n);
607  Map<ArrayXr> mus(m_workspace.data()+2*n, n);
608  Map<ArrayXr> zhat(m_workspace.data()+3*n, n);
609 
610 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
611  std::cout << "computeSVDofM using:\n";
612  std::cout << " z: " << col0.transpose() << "\n";
613  std::cout << " d: " << diag.transpose() << "\n";
614 #endif
615 
616  // Compute singVals, shifts, and mus
617  computeSingVals(col0, diag, perm, singVals, shifts, mus);
618 
619 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
620  std::cout << " j: " << (m_computed.block(firstCol, firstCol, n, n)).jacobiSvd().singularValues().transpose().reverse() << "\n\n";
621  std::cout << " sing-val: " << singVals.transpose() << "\n";
622  std::cout << " mu: " << mus.transpose() << "\n";
623  std::cout << " shift: " << shifts.transpose() << "\n";
624 
625  {
626  std::cout << "\n\n mus: " << mus.head(actual_n).transpose() << "\n\n";
627  std::cout << " check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
628  assert((((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n) >= 0).all());
629  std::cout << " check2 (>0) : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n";
630  assert((((singVals.array()-diag) / singVals.array()).head(actual_n) >= 0).all());
631  }
632 #endif
633 
634 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
635  assert(singVals.allFinite());
636  assert(mus.allFinite());
637  assert(shifts.allFinite());
638 #endif
639 
640  // Compute zhat
641  perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat);
642 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
643  std::cout << " zhat: " << zhat.transpose() << "\n";
644 #endif
645 
646 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
647  assert(zhat.allFinite());
648 #endif
649 
650  computeSingVecs(zhat, diag, perm, singVals, shifts, mus, U, V);
651 
652 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
653  std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n";
654  std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() << "\n";
655 #endif
656 
657 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
658  assert(m_naiveU.allFinite());
659  assert(m_naiveV.allFinite());
660  assert(m_computed.allFinite());
661  assert(U.allFinite());
662  assert(V.allFinite());
663 // assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 100*NumTraits<RealScalar>::epsilon() * n);
664 // assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 100*NumTraits<RealScalar>::epsilon() * n);
665 #endif
666 
667  // Because of deflation, the singular values might not be completely sorted.
668  // Fortunately, reordering them is a O(n) problem
669  for(Index i=0; i<actual_n-1; ++i)
670  {
671  if(singVals(i)>singVals(i+1))
672  {
673  using std::swap;
674  swap(singVals(i),singVals(i+1));
675  U.col(i).swap(U.col(i+1));
676  if(m_compV) V.col(i).swap(V.col(i+1));
677  }
678  }
679 
680 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
681  {
682  bool singular_values_sorted = (((singVals.segment(1,actual_n-1)-singVals.head(actual_n-1))).array() >= 0).all();
683  if(!singular_values_sorted)
684  std::cout << "Singular values are not sorted: " << singVals.segment(1,actual_n).transpose() << "\n";
685  assert(singular_values_sorted);
686  }
687 #endif
688 
689  // Reverse order so that singular values in increased order
690  // Because of deflation, the zeros singular-values are already at the end
691  singVals.head(actual_n).reverseInPlace();
692  U.leftCols(actual_n).rowwise().reverseInPlace();
693  if (m_compV) V.leftCols(actual_n).rowwise().reverseInPlace();
694 
695 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
696  JacobiSVD<MatrixXr> jsvd(m_computed.block(firstCol, firstCol, n, n) );
697  std::cout << " * j: " << jsvd.singularValues().transpose() << "\n\n";
698  std::cout << " * sing-val: " << singVals.transpose() << "\n";
699 // std::cout << " * err: " << ((jsvd.singularValues()-singVals)>1e-13*singVals.norm()).transpose() << "\n";
700 #endif
701 }
702 
703 template <typename MatrixType>
704 typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar mu, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const ArrayRef& diagShifted, RealScalar shift)
705 {
706  Index m = perm.size();
707  RealScalar res = Literal(1);
708  for(Index i=0; i<m; ++i)
709  {
710  Index j = perm(i);
711  // The following expression could be rewritten to involve only a single division,
712  // but this would make the expression more sensitive to overflow.
713  res += (col0(j) / (diagShifted(j) - mu)) * (col0(j) / (diag(j) + shift + mu));
714  }
715  return res;
716 
717 }
718 
719 template <typename MatrixType>
720 void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm,
721  VectorType& singVals, ArrayRef shifts, ArrayRef mus)
722 {
723  using std::abs;
724  using std::swap;
725  using std::sqrt;
726 
727  Index n = col0.size();
728  Index actual_n = n;
729  // Note that here actual_n is computed based on col0(i)==0 instead of diag(i)==0 as above
730  // because 1) we have diag(i)==0 => col0(i)==0 and 2) if col0(i)==0, then diag(i) is already a singular value.
731  while(actual_n>1 && col0(actual_n-1)==Literal(0)) --actual_n;
732 
733  for (Index k = 0; k < n; ++k)
734  {
735  if (col0(k) == Literal(0) || actual_n==1)
736  {
737  // if col0(k) == 0, then entry is deflated, so singular value is on diagonal
738  // if actual_n==1, then the deflated problem is already diagonalized
739  singVals(k) = k==0 ? col0(0) : diag(k);
740  mus(k) = Literal(0);
741  shifts(k) = k==0 ? col0(0) : diag(k);
742  continue;
743  }
744 
745  // otherwise, use secular equation to find singular value
746  RealScalar left = diag(k);
747  RealScalar right; // was: = (k != actual_n-1) ? diag(k+1) : (diag(actual_n-1) + col0.matrix().norm());
748  if(k==actual_n-1)
749  right = (diag(actual_n-1) + col0.matrix().norm());
750  else
751  {
752  // Skip deflated singular values,
753  // recall that at this stage we assume that z[j]!=0 and all entries for which z[j]==0 have been put aside.
754  // This should be equivalent to using perm[]
755  Index l = k+1;
756  while(col0(l)==Literal(0)) { ++l; eigen_internal_assert(l<actual_n); }
757  right = diag(l);
758  }
759 
760  // first decide whether it's closer to the left end or the right end
761  RealScalar mid = left + (right-left) / Literal(2);
762  RealScalar fMid = secularEq(mid, col0, diag, perm, diag, Literal(0));
763 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
764  std::cout << "right-left = " << right-left << "\n";
765 // std::cout << "fMid = " << fMid << " " << secularEq(mid-left, col0, diag, perm, ArrayXr(diag-left), left)
766 // << " " << secularEq(mid-right, col0, diag, perm, ArrayXr(diag-right), right) << "\n";
767  std::cout << " = " << secularEq(left+RealScalar(0.000001)*(right-left), col0, diag, perm, diag, 0)
768  << " " << secularEq(left+RealScalar(0.1) *(right-left), col0, diag, perm, diag, 0)
769  << " " << secularEq(left+RealScalar(0.2) *(right-left), col0, diag, perm, diag, 0)
770  << " " << secularEq(left+RealScalar(0.3) *(right-left), col0, diag, perm, diag, 0)
771  << " " << secularEq(left+RealScalar(0.4) *(right-left), col0, diag, perm, diag, 0)
772  << " " << secularEq(left+RealScalar(0.49) *(right-left), col0, diag, perm, diag, 0)
773  << " " << secularEq(left+RealScalar(0.5) *(right-left), col0, diag, perm, diag, 0)
774  << " " << secularEq(left+RealScalar(0.51) *(right-left), col0, diag, perm, diag, 0)
775  << " " << secularEq(left+RealScalar(0.6) *(right-left), col0, diag, perm, diag, 0)
776  << " " << secularEq(left+RealScalar(0.7) *(right-left), col0, diag, perm, diag, 0)
777  << " " << secularEq(left+RealScalar(0.8) *(right-left), col0, diag, perm, diag, 0)
778  << " " << secularEq(left+RealScalar(0.9) *(right-left), col0, diag, perm, diag, 0)
779  << " " << secularEq(left+RealScalar(0.999999)*(right-left), col0, diag, perm, diag, 0) << "\n";
780 #endif
781  RealScalar shift = (k == actual_n-1 || fMid > Literal(0)) ? left : right;
782 
783  // measure everything relative to shift
784  Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n);
785  diagShifted = diag - shift;
786 
787  // initial guess
788  RealScalar muPrev, muCur;
789  if (shift == left)
790  {
791  muPrev = (right - left) * RealScalar(0.1);
792  if (k == actual_n-1) muCur = right - left;
793  else muCur = (right - left) * RealScalar(0.5);
794  }
795  else
796  {
797  muPrev = -(right - left) * RealScalar(0.1);
798  muCur = -(right - left) * RealScalar(0.5);
799  }
800 
801  RealScalar fPrev = secularEq(muPrev, col0, diag, perm, diagShifted, shift);
802  RealScalar fCur = secularEq(muCur, col0, diag, perm, diagShifted, shift);
803  if (abs(fPrev) < abs(fCur))
804  {
805  swap(fPrev, fCur);
806  swap(muPrev, muCur);
807  }
808 
809  // rational interpolation: fit a function of the form a / mu + b through the two previous
810  // iterates and use its zero to compute the next iterate
811  bool useBisection = fPrev*fCur>Literal(0);
812  while (fCur!=Literal(0) && abs(muCur - muPrev) > Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
813  {
814  ++m_numIters;
815 
816  // Find a and b such that the function f(mu) = a / mu + b matches the current and previous samples.
817  RealScalar a = (fCur - fPrev) / (Literal(1)/muCur - Literal(1)/muPrev);
818  RealScalar b = fCur - a / muCur;
819  // And find mu such that f(mu)==0:
820  RealScalar muZero = -a/b;
821  RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift);
822 
823 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
824  assert((std::isfinite)(fZero));
825 #endif
826 
827  muPrev = muCur;
828  fPrev = fCur;
829  muCur = muZero;
830  fCur = fZero;
831 
832  if (shift == left && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
833  if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
834  if (abs(fCur)>abs(fPrev)) useBisection = true;
835  }
836 
837  // fall back on bisection method if rational interpolation did not work
838  if (useBisection)
839  {
840 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
841  std::cout << "useBisection for k = " << k << ", actual_n = " << actual_n << "\n";
842 #endif
843  RealScalar leftShifted, rightShifted;
844  if (shift == left)
845  {
846  // to avoid overflow, we must have mu > max(real_min, |z(k)|/sqrt(real_max)),
847  // the factor 2 is to be more conservative
848  leftShifted = numext::maxi<RealScalar>( (std::numeric_limits<RealScalar>::min)(), Literal(2) * abs(col0(k)) / sqrt((std::numeric_limits<RealScalar>::max)()) );
849 
850  // check that we did it right:
851  eigen_internal_assert( (numext::isfinite)( (col0(k)/leftShifted)*(col0(k)/(diag(k)+shift+leftShifted)) ) );
852  // I don't understand why the case k==0 would be special there:
853  // if (k == 0) rightShifted = right - left; else
854  rightShifted = (k==actual_n-1) ? right : ((right - left) * RealScalar(0.51)); // theoretically we can take 0.5, but let's be safe
855  }
856  else
857  {
858  leftShifted = -(right - left) * RealScalar(0.51);
859  if(k+1<n)
860  rightShifted = -numext::maxi<RealScalar>( (std::numeric_limits<RealScalar>::min)(), abs(col0(k+1)) / sqrt((std::numeric_limits<RealScalar>::max)()) );
861  else
862  rightShifted = -(std::numeric_limits<RealScalar>::min)();
863  }
864 
865  RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
866 
867 #if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_SANITY_CHECKS
868  RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
869 #endif
870 
871 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
872  if(!(std::isfinite)(fLeft))
873  std::cout << "f(" << leftShifted << ") =" << fLeft << " ; " << left << " " << shift << " " << right << "\n";
874  assert((std::isfinite)(fLeft));
875 
876  if(!(std::isfinite)(fRight))
877  std::cout << "f(" << rightShifted << ") =" << fRight << " ; " << left << " " << shift << " " << right << "\n";
878 // assert((std::isfinite)(fRight));
879 #endif
880 
881 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
882  if(!(fLeft * fRight<0))
883  {
884  std::cout << "f(leftShifted) using leftShifted=" << leftShifted << " ; diagShifted(1:10):" << diagShifted.head(10).transpose() << "\n ; "
885  << "left==shift=" << bool(left==shift) << " ; left-shift = " << (left-shift) << "\n";
886  std::cout << "k=" << k << ", " << fLeft << " * " << fRight << " == " << fLeft * fRight << " ; "
887  << "[" << left << " .. " << right << "] -> [" << leftShifted << " " << rightShifted << "], shift=" << shift << " , f(right)=" << secularEq(0, col0, diag, perm, diagShifted, shift) << " == " << secularEq(right, col0, diag, perm, diag, 0) << "\n";
888  }
889 #endif
890  eigen_internal_assert(fLeft * fRight < Literal(0));
891 
892  while (rightShifted - leftShifted > Literal(2) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted)))
893  {
894  RealScalar midShifted = (leftShifted + rightShifted) / Literal(2);
895  fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
896  eigen_internal_assert((numext::isfinite)(fMid));
897 
898  if (fLeft * fMid < Literal(0))
899  {
900  rightShifted = midShifted;
901  }
902  else
903  {
904  leftShifted = midShifted;
905  fLeft = fMid;
906  }
907  }
908 
909  muCur = (leftShifted + rightShifted) / Literal(2);
910  }
911 
912  singVals[k] = shift + muCur;
913  shifts[k] = shift;
914  mus[k] = muCur;
915 
916 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
917  if(k+1<n)
918  std::cout << "found " << singVals[k] << " == " << shift << " + " << muCur << " from " << diag(k) << " .. " << diag(k+1) << "\n";
919 #endif
920 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
921  assert(k==0 || singVals[k]>=singVals[k-1]);
922  assert(singVals[k]>=diag(k));
923 #endif
924 
925  // perturb singular value slightly if it equals diagonal entry to avoid division by zero later
926  // (deflation is supposed to avoid this from happening)
927  // - this does no seem to be necessary anymore -
928 // if (singVals[k] == left) singVals[k] *= 1 + NumTraits<RealScalar>::epsilon();
929 // if (singVals[k] == right) singVals[k] *= 1 - NumTraits<RealScalar>::epsilon();
930  }
931 }
932 
933 
934 // zhat is perturbation of col0 for which singular vectors can be computed stably (see Section 3.1)
935 template <typename MatrixType>
937  (const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals,
938  const ArrayRef& shifts, const ArrayRef& mus, ArrayRef zhat)
939 {
940  using std::sqrt;
941  Index n = col0.size();
942  Index m = perm.size();
943  if(m==0)
944  {
945  zhat.setZero();
946  return;
947  }
948  Index last = perm(m-1);
949  // The offset permits to skip deflated entries while computing zhat
950  for (Index k = 0; k < n; ++k)
951  {
952  if (col0(k) == Literal(0)) // deflated
953  zhat(k) = Literal(0);
954  else
955  {
956  // see equation (3.6)
957  RealScalar dk = diag(k);
958  RealScalar prod = (singVals(last) + dk) * (mus(last) + (shifts(last) - dk));
959 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
960  if(prod<0) {
961  std::cout << "k = " << k << " ; z(k)=" << col0(k) << ", diag(k)=" << dk << "\n";
962  std::cout << "prod = " << "(" << singVals(last) << " + " << dk << ") * (" << mus(last) << " + (" << shifts(last) << " - " << dk << "))" << "\n";
963  std::cout << " = " << singVals(last) + dk << " * " << mus(last) + (shifts(last) - dk) << "\n";
964  }
965  assert(prod>=0);
966 #endif
967 
968  for(Index l = 0; l<m; ++l)
969  {
970  Index i = perm(l);
971  if(i!=k)
972  {
973 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
974  if(i>=k && (l==0 || l-1>=m))
975  {
976  std::cout << "Error in perturbCol0\n";
977  std::cout << " " << k << "/" << n << " " << l << "/" << m << " " << i << "/" << n << " ; " << col0(k) << " " << diag(k) << " " << "\n";
978  std::cout << " " <<diag(i) << "\n";
979  Index j = (i<k /*|| l==0*/) ? i : perm(l-1);
980  std::cout << " " << "j=" << j << "\n";
981  }
982 #endif
983  Index j = i<k ? i : perm(l-1);
984 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
985  if(!(dk!=Literal(0) || diag(i)!=Literal(0)))
986  {
987  std::cout << "k=" << k << ", i=" << i << ", l=" << l << ", perm.size()=" << perm.size() << "\n";
988  }
989  assert(dk!=Literal(0) || diag(i)!=Literal(0));
990 #endif
991  prod *= ((singVals(j)+dk) / ((diag(i)+dk))) * ((mus(j)+(shifts(j)-dk)) / ((diag(i)-dk)));
992 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
993  assert(prod>=0);
994 #endif
995 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
996  if(i!=k && std::abs(((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) - 1) > 0.9 )
997  std::cout << " " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) << " == (" << (singVals(j)+dk) << " * " << (mus(j)+(shifts(j)-dk))
998  << ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n";
999 #endif
1000  }
1001  }
1002 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1003  std::cout << "zhat(" << k << ") = sqrt( " << prod << ") ; " << (singVals(last) + dk) << " * " << mus(last) + shifts(last) << " - " << dk << "\n";
1004 #endif
1005  RealScalar tmp = sqrt(prod);
1006 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1007  assert((std::isfinite)(tmp));
1008 #endif
1009  zhat(k) = col0(k) > Literal(0) ? tmp : -tmp;
1010  }
1011  }
1012 }
1013 
1014 // compute singular vectors
1015 template <typename MatrixType>
1017  (const ArrayRef& zhat, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals,
1018  const ArrayRef& shifts, const ArrayRef& mus, MatrixXr& U, MatrixXr& V)
1019 {
1020  Index n = zhat.size();
1021  Index m = perm.size();
1022 
1023  for (Index k = 0; k < n; ++k)
1024  {
1025  if (zhat(k) == Literal(0))
1026  {
1027  U.col(k) = VectorType::Unit(n+1, k);
1028  if (m_compV) V.col(k) = VectorType::Unit(n, k);
1029  }
1030  else
1031  {
1032  U.col(k).setZero();
1033  for(Index l=0;l<m;++l)
1034  {
1035  Index i = perm(l);
1036  U(i,k) = zhat(i)/(((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
1037  }
1038  U(n,k) = Literal(0);
1039  U.col(k).normalize();
1040 
1041  if (m_compV)
1042  {
1043  V.col(k).setZero();
1044  for(Index l=1;l<m;++l)
1045  {
1046  Index i = perm(l);
1047  V(i,k) = diag(i) * zhat(i) / (((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
1048  }
1049  V(0,k) = Literal(-1);
1050  V.col(k).normalize();
1051  }
1052  }
1053  }
1054  U.col(n) = VectorType::Unit(n+1, n);
1055 }
1056 
1057 
1058 // page 12_13
1059 // i >= 1, di almost null and zi non null.
1060 // We use a rotation to zero out zi applied to the left of M
1061 template <typename MatrixType>
1062 void BDCSVD<MatrixType>::deflation43(Index firstCol, Index shift, Index i, Index size)
1063 {
1064  using std::abs;
1065  using std::sqrt;
1066  using std::pow;
1067  Index start = firstCol + shift;
1068  RealScalar c = m_computed(start, start);
1069  RealScalar s = m_computed(start+i, start);
1070  RealScalar r = numext::hypot(c,s);
1071  if (r == Literal(0))
1072  {
1073  m_computed(start+i, start+i) = Literal(0);
1074  return;
1075  }
1076  m_computed(start,start) = r;
1077  m_computed(start+i, start) = Literal(0);
1078  m_computed(start+i, start+i) = Literal(0);
1079 
1080  JacobiRotation<RealScalar> J(c/r,-s/r);
1081  if (m_compU) m_naiveU.middleRows(firstCol, size+1).applyOnTheRight(firstCol, firstCol+i, J);
1082  else m_naiveU.applyOnTheRight(firstCol, firstCol+i, J);
1083 }// end deflation 43
1084 
1085 
1086 // page 13
1087 // i,j >= 1, i!=j and |di - dj| < epsilon * norm2(M)
1088 // We apply two rotations to have zj = 0;
1089 // TODO deflation44 is still broken and not properly tested
1090 template <typename MatrixType>
1091 void BDCSVD<MatrixType>::deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size)
1092 {
1093  using std::abs;
1094  using std::sqrt;
1095  using std::conj;
1096  using std::pow;
1097  RealScalar c = m_computed(firstColm+i, firstColm);
1098  RealScalar s = m_computed(firstColm+j, firstColm);
1099  RealScalar r = sqrt(numext::abs2(c) + numext::abs2(s));
1100 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1101  std::cout << "deflation 4.4: " << i << "," << j << " -> " << c << " " << s << " " << r << " ; "
1102  << m_computed(firstColm + i-1, firstColm) << " "
1103  << m_computed(firstColm + i, firstColm) << " "
1104  << m_computed(firstColm + i+1, firstColm) << " "
1105  << m_computed(firstColm + i+2, firstColm) << "\n";
1106  std::cout << m_computed(firstColm + i-1, firstColm + i-1) << " "
1107  << m_computed(firstColm + i, firstColm+i) << " "
1108  << m_computed(firstColm + i+1, firstColm+i+1) << " "
1109  << m_computed(firstColm + i+2, firstColm+i+2) << "\n";
1110 #endif
1111  if (r==Literal(0))
1112  {
1113  m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j);
1114  return;
1115  }
1116  c/=r;
1117  s/=r;
1118  m_computed(firstColm + i, firstColm) = r;
1119  m_computed(firstColm + j, firstColm + j) = m_computed(firstColm + i, firstColm + i);
1120  m_computed(firstColm + j, firstColm) = Literal(0);
1121 
1123  if (m_compU) m_naiveU.middleRows(firstColu, size+1).applyOnTheRight(firstColu + i, firstColu + j, J);
1124  else m_naiveU.applyOnTheRight(firstColu+i, firstColu+j, J);
1125  if (m_compV) m_naiveV.middleRows(firstRowW, size).applyOnTheRight(firstColW + i, firstColW + j, J);
1126 }// end deflation 44
1127 
1128 
1129 // acts on block from (firstCol+shift, firstCol+shift) to (lastCol+shift, lastCol+shift) [inclusive]
1130 template <typename MatrixType>
1131 void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift)
1132 {
1133  using std::sqrt;
1134  using std::abs;
1135  const Index length = lastCol + 1 - firstCol;
1136 
1137  Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1);
1138  Diagonal<MatrixXr> fulldiag(m_computed);
1139  VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length);
1140 
1141  const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
1142  RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
1143  RealScalar epsilon_strict = numext::maxi<RealScalar>(considerZero,NumTraits<RealScalar>::epsilon() * maxDiag);
1144  RealScalar epsilon_coarse = Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag);
1145 
1146 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1147  assert(m_naiveU.allFinite());
1148  assert(m_naiveV.allFinite());
1149  assert(m_computed.allFinite());
1150 #endif
1151 
1152 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1153  std::cout << "\ndeflate:" << diag.head(k+1).transpose() << " | " << diag.segment(k+1,length-k-1).transpose() << "\n";
1154 #endif
1155 
1156  //condition 4.1
1157  if (diag(0) < epsilon_coarse)
1158  {
1159 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1160  std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon_coarse << "\n";
1161 #endif
1162  diag(0) = epsilon_coarse;
1163  }
1164 
1165  //condition 4.2
1166  for (Index i=1;i<length;++i)
1167  if (abs(col0(i)) < epsilon_strict)
1168  {
1169 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1170  std::cout << "deflation 4.2, set z(" << i << ") to zero because " << abs(col0(i)) << " < " << epsilon_strict << " (diag(" << i << ")=" << diag(i) << ")\n";
1171 #endif
1172  col0(i) = Literal(0);
1173  }
1174 
1175  //condition 4.3
1176  for (Index i=1;i<length; i++)
1177  if (diag(i) < epsilon_coarse)
1178  {
1179 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1180  std::cout << "deflation 4.3, cancel z(" << i << ")=" << col0(i) << " because diag(" << i << ")=" << diag(i) << " < " << epsilon_coarse << "\n";
1181 #endif
1182  deflation43(firstCol, shift, i, length);
1183  }
1184 
1185 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1186  assert(m_naiveU.allFinite());
1187  assert(m_naiveV.allFinite());
1188  assert(m_computed.allFinite());
1189 #endif
1190 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1191  std::cout << "to be sorted: " << diag.transpose() << "\n\n";
1192  std::cout << " : " << col0.transpose() << "\n\n";
1193 #endif
1194  {
1195  // Check for total deflation
1196  // If we have a total deflation, then we have to consider col0(0)==diag(0) as a singular value during sorting
1197  bool total_deflation = (col0.tail(length-1).array()<considerZero).all();
1198 
1199  // Sort the diagonal entries, since diag(1:k-1) and diag(k:length) are already sorted, let's do a sorted merge.
1200  // First, compute the respective permutation.
1201  Index *permutation = m_workspaceI.data();
1202  {
1203  permutation[0] = 0;
1204  Index p = 1;
1205 
1206  // Move deflated diagonal entries at the end.
1207  for(Index i=1; i<length; ++i)
1208  if(abs(diag(i))<considerZero)
1209  permutation[p++] = i;
1210 
1211  Index i=1, j=k+1;
1212  for( ; p < length; ++p)
1213  {
1214  if (i > k) permutation[p] = j++;
1215  else if (j >= length) permutation[p] = i++;
1216  else if (diag(i) < diag(j)) permutation[p] = j++;
1217  else permutation[p] = i++;
1218  }
1219  }
1220 
1221  // If we have a total deflation, then we have to insert diag(0) at the right place
1222  if(total_deflation)
1223  {
1224  for(Index i=1; i<length; ++i)
1225  {
1226  Index pi = permutation[i];
1227  if(abs(diag(pi))<considerZero || diag(0)<diag(pi))
1228  permutation[i-1] = permutation[i];
1229  else
1230  {
1231  permutation[i-1] = 0;
1232  break;
1233  }
1234  }
1235  }
1236 
1237  // Current index of each col, and current column of each index
1238  Index *realInd = m_workspaceI.data()+length;
1239  Index *realCol = m_workspaceI.data()+2*length;
1240 
1241  for(int pos = 0; pos< length; pos++)
1242  {
1243  realCol[pos] = pos;
1244  realInd[pos] = pos;
1245  }
1246 
1247  for(Index i = total_deflation?0:1; i < length; i++)
1248  {
1249  const Index pi = permutation[length - (total_deflation ? i+1 : i)];
1250  const Index J = realCol[pi];
1251 
1252  using std::swap;
1253  // swap diagonal and first column entries:
1254  swap(diag(i), diag(J));
1255  if(i!=0 && J!=0) swap(col0(i), col0(J));
1256 
1257  // change columns
1258  if (m_compU) m_naiveU.col(firstCol+i).segment(firstCol, length + 1).swap(m_naiveU.col(firstCol+J).segment(firstCol, length + 1));
1259  else m_naiveU.col(firstCol+i).segment(0, 2) .swap(m_naiveU.col(firstCol+J).segment(0, 2));
1260  if (m_compV) m_naiveV.col(firstColW + i).segment(firstRowW, length).swap(m_naiveV.col(firstColW + J).segment(firstRowW, length));
1261 
1262  //update real pos
1263  const Index realI = realInd[i];
1264  realCol[realI] = J;
1265  realCol[pi] = i;
1266  realInd[J] = realI;
1267  realInd[i] = pi;
1268  }
1269  }
1270 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1271  std::cout << "sorted: " << diag.transpose().format(bdcsvdfmt) << "\n";
1272  std::cout << " : " << col0.transpose() << "\n\n";
1273 #endif
1274 
1275  //condition 4.4
1276  {
1277  Index i = length-1;
1278  while(i>0 && (abs(diag(i))<considerZero || abs(col0(i))<considerZero)) --i;
1279  for(; i>1;--i)
1280  if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag )
1281  {
1282 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1283  std::cout << "deflation 4.4 with i = " << i << " because " << diag(i) << " - " << diag(i-1) << " == " << (diag(i) - diag(i-1)) << " < " << NumTraits<RealScalar>::epsilon()*/*diag(i)*/maxDiag << "\n";
1284 #endif
1285  eigen_internal_assert(abs(diag(i) - diag(i-1))<epsilon_coarse && " diagonal entries are not properly sorted");
1286  deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length);
1287  }
1288  }
1289 
1290 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1291  for(Index j=2;j<length;++j)
1292  assert(diag(j-1)<=diag(j) || abs(diag(j))<considerZero);
1293 #endif
1294 
1295 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1296  assert(m_naiveU.allFinite());
1297  assert(m_naiveV.allFinite());
1298  assert(m_computed.allFinite());
1299 #endif
1300 }//end deflation
1301 
1302 #ifndef EIGEN_CUDACC
1303 
1309 template<typename Derived>
1311 MatrixBase<Derived>::bdcSvd(unsigned int computationOptions) const
1312 {
1313  return BDCSVD<PlainObject>(*this, computationOptions);
1314 }
1315 #endif
1316 
1317 } // end namespace Eigen
1318 
1319 #endif
BDCSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: BDCSVD.h:167
static const Eigen::internal::all_t all
Definition: IndexedViewHelper.h:181
Definition: Constants.h:387
BDCSVD< PlainObject > bdcSvd(unsigned int computationOptions=0) const
Definition: BDCSVD.h:1311
const AdjointReturnType adjoint() const
Definition: Transpose.h:210
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
const MatrixUType & matrixU() const
Definition: SVDBase.h:83
Namespace containing all symbols from the Eigen library.
Definition: Core:324
Rotation given by a cosine-sine pair.
Definition: ForwardDeclarations.h:264
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:184
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:279
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515
Expression of a fixed-size or dynamic-size sub-vector.
Definition: ForwardDeclarations.h:88
BDCSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: BDCSVD.h:123
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Definition: Constants.h:239
BDCSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: BDCSVD.h:243
BDCSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: BDCSVD.h:139
class Bidiagonal Divide and Conquer SVD
Definition: ForwardDeclarations.h:260
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:192
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
Index nonzeroSingularValues() const
Definition: SVDBase.h:118
BDCSVD()
Default Constructor.
Definition: BDCSVD.h:113
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Definition: ForwardDeclarations.h:259
const SingularValuesType & singularValues() const
Definition: SVDBase.h:111
const MatrixVType & matrixV() const
Definition: SVDBase.h:99
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:63
const int Dynamic
Definition: Constants.h:21
Definition: Constants.h:391
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
static const Symbolic::SymbolExpr< internal::symbolic_last_tag > last
Definition: IndexedViewHelper.h:45