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> |