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