This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 319
Collapse All | Expand All

(-)a/Eigen/src/Cholesky/LDLT.h (-67 / +44 lines)
Lines 90-126 template<typename _MatrixType, int _UpLo Link Here
90
      * Like the default constructor but with preallocation of the internal data
90
      * Like the default constructor but with preallocation of the internal data
91
      * according to the specified problem \a size.
91
      * according to the specified problem \a size.
92
      * \sa LDLT()
92
      * \sa LDLT()
93
      */
93
      */
94
    LDLT(Index size)
94
    LDLT(Index size)
95
      : m_matrix(size, size),
95
      : m_matrix(size, size),
96
        m_transpositions(size),
96
        m_transpositions(size),
97
        m_temporary(size),
97
        m_temporary(size),
98
        m_temporary2(size),
99
        m_isInitialized(false)
98
        m_isInitialized(false)
100
    {}
99
    {}
101
100
102
    /** \brief Constructor with decomposition
101
    /** \brief Constructor with decomposition
103
      *
102
      *
104
      * This calculates the decomposition for the input \a matrix.
103
      * This calculates the decomposition for the input \a matrix.
105
      * \sa LDLT(Index size)
104
      * \sa LDLT(Index size)
106
      */
105
      */
107
    LDLT(const MatrixType& matrix)
106
    LDLT(const MatrixType& matrix)
108
      : m_matrix(matrix.rows(), matrix.cols()),
107
      : m_matrix(matrix.rows(), matrix.cols()),
109
        m_transpositions(matrix.rows()),
108
        m_transpositions(matrix.rows()),
110
        m_temporary(matrix.rows()),
109
        m_temporary(matrix.rows()),
111
        m_temporary2(matrix.rows()),
112
        m_isInitialized(false)
110
        m_isInitialized(false)
113
    {
111
    {
114
      compute(matrix);
112
      compute(matrix);
115
    }
113
    }
116
114
117
    /** Clear any existing decomposition
115
    /** Clear any existing decomposition
118
     * \sa update(W,sigma)
116
     * \sa rankUpdate(w,sigma)
119
     */
117
     */
120
    void clear()
118
    void clear()
121
    {
119
    {
122
      m_isInitialized = false;
120
      m_isInitialized = false;
123
    }
121
    }
124
122
125
    /** \returns a view of the upper triangular matrix U */
123
    /** \returns a view of the upper triangular matrix U */
126
    inline typename Traits::MatrixU matrixU() const
124
    inline typename Traits::MatrixU matrixU() const
Lines 197-215 template<typename _MatrixType, int _UpLo Link Here
197
    }
195
    }
198
    #endif
196
    #endif
199
197
200
    template<typename Derived>
198
    template<typename Derived>
201
    bool solveInPlace(MatrixBase<Derived> &bAndX) const;
199
    bool solveInPlace(MatrixBase<Derived> &bAndX) const;
202
200
203
    LDLT& compute(const MatrixType& matrix);
201
    LDLT& compute(const MatrixType& matrix);
204
202
205
    // Note W is modified, this is just a hack to allow, e.g., W.col(i)
206
    template <typename Derived>
203
    template <typename Derived>
207
    LDLT& update(MatrixBase<Derived>& W,int sigma=1);
204
    LDLT& rankUpdate(const MatrixBase<Derived>& w,RealScalar alpha=1);
208
205
209
    /** \returns the internal LDLT decomposition matrix
206
    /** \returns the internal LDLT decomposition matrix
210
      *
207
      *
211
      * TODO: document the storage layout
208
      * TODO: document the storage layout
212
      */
209
      */
213
    inline const MatrixType& matrixLDLT() const
210
    inline const MatrixType& matrixLDLT() const
214
    {
211
    {
215
      eigen_assert(m_isInitialized && "LDLT is not initialized.");
212
      eigen_assert(m_isInitialized && "LDLT is not initialized.");
Lines 226-242 template<typename _MatrixType, int _UpLo Link Here
226
    /** \internal
223
    /** \internal
227
      * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U.
224
      * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U.
228
      * The strict upper part is used during the decomposition, the strict lower
225
      * The strict upper part is used during the decomposition, the strict lower
229
      * part correspond to the coefficients of L (its diagonal is equal to 1 and
226
      * part correspond to the coefficients of L (its diagonal is equal to 1 and
230
      * is not stored), and the diagonal entries correspond to D.
227
      * is not stored), and the diagonal entries correspond to D.
231
      */
228
      */
232
    MatrixType m_matrix;
229
    MatrixType m_matrix;
233
    TranspositionType m_transpositions;
230
    TranspositionType m_transpositions;
234
    TmpMatrixType m_temporary, m_temporary2;
231
    TmpMatrixType m_temporary;
235
    int m_sign;
232
    int m_sign;
236
    bool m_isInitialized;
233
    bool m_isInitialized;
237
};
234
};
238
235
239
namespace internal {
236
namespace internal {
240
237
241
template<int UpLo> struct ldlt_inplace;
238
template<int UpLo> struct ldlt_inplace;
242
239
Lines 323-429 template<> struct ldlt_inplace<Lower> Link Here
323
      }
320
      }
324
      if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
321
      if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
325
        A21 /= mat.coeffRef(k,k);
322
        A21 /= mat.coeffRef(k,k);
326
    }
323
    }
327
324
328
    return true;
325
    return true;
329
  }
326
  }
330
327
331
  // Note: if the original matrix is not of full rank, you can only do
328
  // Reference for the algorithm: Davis and Hager, "Multiple Rank
332
  // rank-1 updates. See update.
329
  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
333
  // Reference for the algorithm: Davis and Hager, "Multiple Rank Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
330
  // Trivial rearrangements of their computations (Timothy E. Holy)
334
  // Trivial rearrangements of their computations (Timothy E. Holy) allow their algorithm to work for rank-1 updates even if the original matrix is not of full rank. 
331
  // allow their algorithm to work for rank-1 updates even if the
335
  template<typename MatrixType, typename Workspace, typename WDerived>
332
  // original matrix is not of full rank.
336
  static bool update_work(MatrixType& mat, Workspace& alpha, Workspace& gamma, MatrixBase<WDerived> const& W, int sigma=1)
333
  // Here only rank-1 updates are implemented, to reduce the
334
  // requirement for intermediate storage and improve accuracy
335
  template<typename MatrixType, typename WDerived>
336
  static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, typename MatrixType::RealScalar sigma=1)
337
  {
337
  {
338
    typedef typename MatrixType::Scalar Scalar;
338
    typedef typename MatrixType::Scalar Scalar;
339
    typedef typename MatrixType::RealScalar RealScalar;
339
    typedef typename MatrixType::RealScalar RealScalar;
340
    typedef typename MatrixType::Index Index;
340
    typedef typename MatrixType::Index Index;
341
    const Index size = mat.rows();
341
    const Index size = mat.rows();
342
    const Index r = W.cols();  // rank of update
342
343
    eigen_assert(mat.cols() == size && w.size()==size);
343
344
344
    // Prepare the update
345
    // Prepare the update
345
    RealScalar alphabar,temp,dtemp,wtemp,denom;
346
    RealScalar alpha,alphabar,temp,dtemp,gammatmp;
346
    for (Index i = 0; i < r; i++)
347
    Scalar wtemp,gamma;
347
      alpha.coeffRef(i) = 1;
348
349
    alpha = 1;
348
350
349
    // Apply the update
351
    // Apply the update
350
    for (Index j = 0; j < size; j++)
352
    for (Index j = 0; j < size; j++)
351
    {
353
    {
352
      // Check for termination due to an original decomposition of low-rank
354
      // Check for termination due to an original decomposition of low-rank
353
      if (!isfinite(alpha.coeff(0)))
355
      if (!std::isfinite(alpha))
354
      {
355
	if (j < r)
356
	  return false;  // failed, too many updates given rank of starting matrix
357
	break;
356
	break;
358
      }
359
357
360
      // Update the diagonal terms
358
      // Update the diagonal terms
361
      dtemp = real(mat.diagonal().coeff(j));
359
      dtemp = real(mat.diagonal().coeff(j));
362
      for (Index i = 0; i < r; i++)
360
      wtemp = w.coeff(j);
363
      {
361
      temp = sigma*real(wtemp*conj(wtemp));
364
	wtemp = W.coeff(j,i);
362
      alphabar = alpha + temp/dtemp;
365
	temp = sigma*wtemp*conj(wtemp);
363
      gammatmp = dtemp*alpha + temp;
366
	alphabar = alpha.coeff(i) + temp/dtemp;
364
      if (gammatmp != 0)
367
	denom = dtemp*alpha.coeff(i) + temp;
365
	gamma = conj(wtemp)/gammatmp;  // FIXME: guessing on conj here
368
	if (denom != 0)
366
      else
369
	  gamma.coeffRef(i) = conj(wtemp)/denom;  // FIXME: guessing on conj here
367
	gamma = 0;
370
	else
368
      dtemp += temp/alpha;
371
	  gamma.coeffRef(i) = 0;
369
      alpha = alphabar;
372
	dtemp += temp/alpha.coeff(i);
373
	alpha.coeffRef(i) = alphabar;
374
      }
375
      mat.diagonal().coeffRef(j) = dtemp;
370
      mat.diagonal().coeffRef(j) = dtemp;
371
376
      // Update the terms of L
372
      // Update the terms of L
377
      for (Index i = 0; i < r; i++)
373
      w.tail(size-j-1) -= wtemp*mat.col(j).tail(size-j-1);
378
      {
374
      mat.col(j).tail(size-j-1) += (sigma*gamma)*w.tail(size-j-1);
379
	wtemp = W.coeff(j,i);
380
	const_cast<MatrixBase<WDerived>& >(W).col(i).tail(size-j-1) -= wtemp*mat.col(j).tail(size-j-1);
381
	mat.col(j).tail(size-j-1) += sigma*gamma.coeff(i)*W.col(i).tail(size-j-1);
382
      }
383
    }
375
    }
384
    return true;
376
    return true;
385
  }
377
  }
386
378
387
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
379
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
388
  static bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& alpha, Workspace& gamma, WType& W, int sigma=1)
380
  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, typename MatrixType::RealScalar sigma=1)
389
  {
381
  {
390
    typedef typename MatrixType::Index Index;
382
    // Apply the permutation to the input w
391
    const Index r = W.cols();  // rank of update
383
    tmp = transpositions * w;
392
384
393
    // Apply the permutation to the input W
385
    return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
394
    W = transpositions * W;
395
396
    // Check to see whether the existing decomposition is
397
    // rank-deficient; if so, do only rank-1 updates
398
    if ((mat.diagonal().array() == 0).any())
399
    {
400
      for (Index i = 0; i < r; i++)
401
	ldlt_inplace<Lower>::update_work(mat,alpha,gamma,W.col(i),sigma);
402
      return true;
403
    } else
404
      return ldlt_inplace<Lower>::update_work(mat,alpha,gamma,W,sigma);
405
  }
386
  }
406
};
387
};
407
388
408
template<> struct ldlt_inplace<Upper>
389
template<> struct ldlt_inplace<Upper>
409
{
390
{
410
  template<typename MatrixType, typename TranspositionType, typename Workspace>
391
  template<typename MatrixType, typename TranspositionType, typename Workspace>
411
  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
392
  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
412
  {
393
  {
413
    Transpose<MatrixType> matt(mat);
394
    Transpose<MatrixType> matt(mat);
414
    return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
395
    return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
415
  }
396
  }
416
397
417
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
398
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
418
  static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& alpha, Workspace& gamma, WType& W, int sigma=1)
399
  static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, typename MatrixType::RealScalar sigma=1)
419
  {
400
  {
420
    Transpose<MatrixType> matt(mat);
401
    Transpose<MatrixType> matt(mat);
421
    return ldlt_inplace<Upper>::update(matt, transpositions, alpha, gamma, W, sigma);
402
    return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w, sigma);
422
  }
403
  }
423
};
404
};
424
405
425
template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
406
template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
426
{
407
{
427
  typedef TriangularView<MatrixType, UnitLower> MatrixL;
408
  typedef TriangularView<MatrixType, UnitLower> MatrixL;
428
  typedef TriangularView<typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
409
  typedef TriangularView<typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
429
  inline static MatrixL getL(const MatrixType& m) { return m; }
410
  inline static MatrixL getL(const MatrixType& m) { return m; }
Lines 448-499 LDLT<MatrixType,_UpLo>& LDLT<MatrixType, Link Here
448
  eigen_assert(a.rows()==a.cols());
429
  eigen_assert(a.rows()==a.cols());
449
  const Index size = a.rows();
430
  const Index size = a.rows();
450
431
451
  m_matrix = a;
432
  m_matrix = a;
452
433
453
  m_transpositions.resize(size);
434
  m_transpositions.resize(size);
454
  m_isInitialized = false;
435
  m_isInitialized = false;
455
  m_temporary.resize(size);
436
  m_temporary.resize(size);
456
  m_temporary2.resize(size);
457
437
458
  internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
438
  internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
459
439
460
  m_isInitialized = true;
440
  m_isInitialized = true;
461
  return *this;
441
  return *this;
462
}
442
}
463
443
464
/** Update the LDLT decomposition:  given A = L D L^T, efficiently compute the decomposition of A + sigma W W^T.  Most commonly, W will be a single column vector, corresponding to a rank-1 update.  However, in general W may be a matrix with multiple columns, although it cannot have more columns than rows.
444
/** Update the LDLT decomposition:  given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T.
465
 * \param W The matrix whose columns are to be incorporated into the decomposition.
445
 * \param w a vector to be incorporated into the decomposition.
466
 * \param sigma An integer, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1.
446
 * \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1.
467
 * \note This function overwrites W.
468
 * \sa clear()
447
 * \sa clear()
469
  */
448
  */
470
template<typename MatrixType, int _UpLo>
449
template<typename MatrixType, int _UpLo>
471
template<typename Derived>
450
template<typename Derived>
472
LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::update(MatrixBase<Derived> & W,int sigma)
451
LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w,typename NumTraits<typename MatrixType::Scalar>::Real sigma)
473
{
452
{
474
  const Index size = W.rows();
453
  const Index size = w.rows();
475
  if (m_isInitialized)
454
  if (m_isInitialized)
476
    eigen_assert(m_matrix.rows()==size);
455
    eigen_assert(m_matrix.rows()==size);
477
  else
456
  else
478
  {    
457
  {    
479
    m_matrix.resize(size,size);
458
    m_matrix.resize(size,size);
480
    m_matrix.setZero();
459
    m_matrix.setZero();
481
    m_transpositions.resize(size);
460
    m_transpositions.resize(size);
482
    for (Index i = 0; i < size; i++)
461
    for (Index i = 0; i < size; i++)
483
      m_transpositions.coeffRef(i) = i;
462
      m_transpositions.coeffRef(i) = i;
484
    m_temporary.resize(size);
463
    m_temporary.resize(size);
485
    m_temporary2.resize(size);
486
    m_sign = sigma;
464
    m_sign = sigma;
487
    m_isInitialized = true;
465
    m_isInitialized = true;
488
  }
466
  }
489
  eigen_assert(W.cols() <= size);  // solely so that temporaries are of adequate size
490
467
491
  internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, m_temporary2, W, sigma);
468
  internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
492
469
493
  return *this;
470
  return *this;
494
}
471
}
495
472
496
namespace internal {
473
namespace internal {
497
template<typename _MatrixType, int _UpLo, typename Rhs>
474
template<typename _MatrixType, int _UpLo, typename Rhs>
498
struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
475
struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
499
  : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
476
  : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
(-)a/test/cholesky.cpp (-1 / +18 lines)
Lines 182-198 template<typename MatrixType> void chole Link Here
182
  VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
182
  VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
183
  m2 = m1;
183
  m2 = m1;
184
  m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB);
184
  m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB);
185
  VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
185
  VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
186
  m2 = m1;
186
  m2 = m1;
187
  m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
187
  m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
188
  VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
188
  VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
189
189
190
  // Cholesky update/downdate
190
  // LLT update/downdate
191
  {
191
  {
192
    MatrixType symmLo = symm.template triangularView<Lower>();
192
    MatrixType symmLo = symm.template triangularView<Lower>();
193
    MatrixType symmUp = symm.template triangularView<Upper>();
193
    MatrixType symmUp = symm.template triangularView<Upper>();
194
    
194
    
195
    VectorType vec = VectorType::Random(rows);
195
    VectorType vec = VectorType::Random(rows);
196
196
197
    MatrixType symmCpy = symm + vec * vec.adjoint();
197
    MatrixType symmCpy = symm + vec * vec.adjoint();
198
198
Lines 200-215 template<typename MatrixType> void chole Link Here
200
    chollo.rankUpdate(vec);
200
    chollo.rankUpdate(vec);
201
    VERIFY_IS_APPROX(symmCpy,  chollo.reconstructedMatrix());
201
    VERIFY_IS_APPROX(symmCpy,  chollo.reconstructedMatrix());
202
202
203
    LLT<MatrixType,Upper> cholup(symmUp);
203
    LLT<MatrixType,Upper> cholup(symmUp);
204
    cholup.rankUpdate(vec);
204
    cholup.rankUpdate(vec);
205
    VERIFY_IS_APPROX(symmCpy,  cholup.reconstructedMatrix());
205
    VERIFY_IS_APPROX(symmCpy,  cholup.reconstructedMatrix());
206
  }
206
  }
207
  
207
  
208
  // LDLT update/downdate
209
  {
210
    MatrixType symmLo = symm.template triangularView<Lower>();
211
    MatrixType symmUp = symm.template triangularView<Upper>();
212
    
213
    VectorType vec = VectorType::Random(rows);
214
215
    MatrixType symmCpy = symm + vec * vec.adjoint();
216
217
    LDLT<MatrixType,Lower> chollo(symmLo);
218
    chollo.rankUpdate(vec);
219
    VERIFY_IS_APPROX(symmCpy,  chollo.reconstructedMatrix());
220
221
    LDLT<MatrixType,Upper> cholup(symmUp);
222
    cholup.rankUpdate(vec);
223
    VERIFY_IS_APPROX(symmCpy,  cholup.reconstructedMatrix());
224
  }
208
}
225
}
209
226
210
template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
227
template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
211
{
228
{
212
  // classic test
229
  // classic test
213
  cholesky(m);
230
  cholesky(m);
214
231
215
  // test mixing real/scalar types
232
  // test mixing real/scalar types

Return to bug 319