Attachment #341 for bug #608

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

View | Details | Raw Unified | Return to bug 608
Collapse All | Expand All




    typedef typename MatrixType::Index Index;
    eigen_assert(mat.rows()==mat.cols());
    const Index size = mat.rows();

    if (size <= 1)
    {
      transpositions.setIdentity();
      if(sign)
        *sign = math::real(mat.coeff(0,0))>0 ? 1:-1;
      return true;
    }

    RealScalar cutoff(0), biggest_in_corner;

    for (Index k = 0; k < size; ++k)
    {
      // Find largest diagonal element

      index_of_biggest_in_corner += k;

      if(k == 0)
      {
        // The biggest overall is the point of reference to which further diagonals
        // are compared; if any diagonal is negligible compared
        // to the largest overall, the algorithm bails.
        cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);

        if(sign)
          *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
      }
      else if(sign)
      {
        // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
        int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
        if(newSign != *sign)
          *sign = 0;
      }

      // Finish early if the matrix is not full rank.
      if(biggest_in_corner < cutoff)
      {
        for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
        if(sign) *sign = 0;
        break;
      }

      transpositions.coeffRef(k) = index_of_biggest_in_corner;
      if(k != index_of_biggest_in_corner)
      {
        // apply the transposition while taking care to consider only
        // the lower triangular part
        Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
        mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
        mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
        std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
        for(int i=k+1;i<index_of_biggest_in_corner;++i)
        {
          Scalar tmp = mat.coeffRef(i,k);
          mat.coeffRef(i,k) = math::conj(mat.coeffRef(index_of_biggest_in_corner,i));
          mat.coeffRef(index_of_biggest_in_corner,i) = math::conj(tmp);
        }
        if(NumTraits<Scalar>::IsComplex)
          mat.coeffRef(index_of_biggest_in_corner,k) = math::conj(mat.coeff(index_of_biggest_in_corner,k));
      }

      // partition the matrix:
      //       A00 |  -  |  -
      // lu  = A10 | A11 |  -
      //       A20 | A21 | A22
      Index rs = size - k - 1;
      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);

      {
        temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
        mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
        if(rs>0)
          A21.noalias() -= A20 * temp.head(k);
      }
      if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
        A21 /= mat.coeffRef(k,k);
      
      if(sign)
      {
        // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
        int newSign = math::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
        if(k == 0)
          *sign = newSign;
        else if(*sign != newSign)
          *sign = 0;
      }
    }

    return true;
  }

  // Reference for the algorithm: Davis and Hager, "Multiple Rank
  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
  // Trivial rearrangements of their computations (Timothy E. Holy)

    // Apply the update
    for (Index j = 0; j < size; j++)
    {
      // Check for termination due to an original decomposition of low-rank
      if (!(isfinite)(alpha))
        break;

      // Update the diagonal terms
      RealScalar dj = math::real(mat.coeff(j,j));
      Scalar wj = w.coeff(j);
      RealScalar swj2 = sigma*abs2(wj);
      RealScalar gamma = dj*alpha + swj2;

      mat.coeffRef(j,j) += swj2/alpha;
      alpha += swj2/dj;


      // Update the terms of L
      Index rs = size-j-1;
      w.tail(rs) -= wj * mat.col(j).tail(rs);
      if(gamma != 0)
        mat.col(j).tail(rs) += (sigma*math::conj(wj)/gamma)*w.tail(rs);
    }
    return true;
  }

  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
  {
    // Apply the permutation to the input w




    }
  }
  else
  {
    temp = vec;
    RealScalar beta = 1;
    for(Index j=0; j<n; ++j)
    {
      RealScalar Ljj = math::real(mat.coeff(j,j));
      RealScalar dj = abs2(Ljj);
      Scalar wj = temp.coeff(j);
      RealScalar swj2 = sigma*abs2(wj);
      RealScalar gamma = dj*beta + swj2;

      RealScalar x = dj + swj2/beta;
      if (x<=RealScalar(0))
        return j;

      beta += swj2/dj;

      // Update the terms of L
      Index rs = n-j-1;
      if(rs)
      {
        temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
        if(gamma != 0)
          mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*math::conj(wj)/gamma)*temp.tail(rs);
      }
    }
  }
  return -1;
}

template<typename Scalar> struct llt_inplace<Scalar, Lower>
{

    for(Index k = 0; k < size; ++k)
    {
      Index rs = size-k-1; // remaining size

      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
      Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
      Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);

      RealScalar x = math::real(mat.coeff(k,k));
      if (k>0) x -= A10.squaredNorm();
      if (x<=RealScalar(0))
        return k;
      mat.coeffRef(k,k) = x = sqrt(x);
      if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
      if (rs>0) A21 *= RealScalar(1)/x;
    }
    return -1;





// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
// This test checks that LDLT reports correctly that matrix is indefinite. 
// See http://forum.kde.org/viewtopic.php?f=74&t=106942
template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
{
  eigen_assert(m.rows() == 2 && m.cols() == 2);
  MatrixType mat;
  {
    mat << 1, 0, 0, -1;
    LDLT<MatrixType> ldlt(mat);
    VERIFY(!ldlt.isNegative());
    VERIFY(!ldlt.isPositive());
  }
  {
    mat << 1, 2, 2, 1;
    LDLT<MatrixType> ldlt(mat);
    VERIFY(!ldlt.isNegative());
    VERIFY(!ldlt.isPositive());
  }
}

template<typename MatrixType> void cholesky_verify_assert()
{
  MatrixType tmp;

  LLT<MatrixType> llt;
  VERIFY_RAISES_ASSERT(llt.matrixL())

Return to bug 608

Lines 254-270 template<> struct ldlt_inplace<Lower> Link Here

(-)a/Eigen/src/Cholesky/LDLT.h (-16 / +17 lines)
254	typedef typename MatrixType::Index Index;	254	typedef typename MatrixType::Index Index;
255	eigen_assert(mat.rows()==mat.cols());	255	eigen_assert(mat.rows()==mat.cols());
256	const Index size = mat.rows();	256	const Index size = mat.rows();
257		257
258	if (size <= 1)	258	if (size <= 1)
259	{	259	{
260	transpositions.setIdentity();	260	transpositions.setIdentity();
261	if(sign)	261	if(sign)
262	*sign = real(mat.coeff(0,0))>0 ? 1:-1;	262	*sign = math::real(mat.coeff(0,0))>0 ? 1:-1;
263	return true;	263	return true;
264	}	264	}
265		265
266	RealScalar cutoff(0), biggest_in_corner;	266	RealScalar cutoff(0), biggest_in_corner;
267		267
268	for (Index k = 0; k < size; ++k)	268	for (Index k = 0; k < size; ++k)
269	{	269	{
270	// Find largest diagonal element	270	// Find largest diagonal element
Lines 273-324 template<> struct ldlt_inplace<Lower> Link Here
273	index_of_biggest_in_corner += k;	273	index_of_biggest_in_corner += k;
274		274
275	if(k == 0)	275	if(k == 0)
276	{	276	{
277	// The biggest overall is the point of reference to which further diagonals	277	// The biggest overall is the point of reference to which further diagonals
278	// are compared; if any diagonal is negligible compared	278	// are compared; if any diagonal is negligible compared
279	// to the largest overall, the algorithm bails.	279	// to the largest overall, the algorithm bails.
280	cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);	280	cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
281
282	if(sign)
283	*sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
284	}
285	else if(sign)
286	{
287	// LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
288	int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
289	if(newSign != *sign)
290	*sign = 0;
291	}	281	}
292		282
293	// Finish early if the matrix is not full rank.	283	// Finish early if the matrix is not full rank.
294	if(biggest_in_corner < cutoff)	284	if(biggest_in_corner < cutoff)
295	{	285	{
296	for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;	286	for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
		287	if(sign) *sign = 0;
297	break;	288	break;
298	}	289	}
299		290
300	transpositions.coeffRef(k) = index_of_biggest_in_corner;	291	transpositions.coeffRef(k) = index_of_biggest_in_corner;
301	if(k != index_of_biggest_in_corner)	292	if(k != index_of_biggest_in_corner)
302	{	293	{
303	// apply the transposition while taking care to consider only	294	// apply the transposition while taking care to consider only
304	// the lower triangular part	295	// the lower triangular part
305	Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element	296	Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
306	mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));	297	mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
307	mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));	298	mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
308	std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));	299	std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
309	for(int i=k+1;i<index_of_biggest_in_corner;++i)	300	for(int i=k+1;i<index_of_biggest_in_corner;++i)
310	{	301	{
311	Scalar tmp = mat.coeffRef(i,k);	302	Scalar tmp = mat.coeffRef(i,k);
312	mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));	303	mat.coeffRef(i,k) = math::conj(mat.coeffRef(index_of_biggest_in_corner,i));
313	mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);	304	mat.coeffRef(index_of_biggest_in_corner,i) = math::conj(tmp);
314	}	305	}
315	if(NumTraits<Scalar>::IsComplex)	306	if(NumTraits<Scalar>::IsComplex)
316	mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));	307	mat.coeffRef(index_of_biggest_in_corner,k) = math::conj(mat.coeff(index_of_biggest_in_corner,k));
317	}	308	}
318		309
319	// partition the matrix:	310	// partition the matrix:
320	// A00 \| - \| -	311	// A00 \| - \| -
321	// lu = A10 \| A11 \| -	312	// lu = A10 \| A11 \| -
322	// A20 \| A21 \| A22	313	// A20 \| A21 \| A22
323	Index rs = size - k - 1;	314	Index rs = size - k - 1;
324	Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);	315	Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
Lines 329-344 template<> struct ldlt_inplace<Lower> Link Here
329	{	320	{
330	temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();	321	temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
331	mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();	322	mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
332	if(rs>0)	323	if(rs>0)
333	A21.noalias() -= A20 * temp.head(k);	324	A21.noalias() -= A20 * temp.head(k);
334	}	325	}
335	if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))	326	if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
336	A21 /= mat.coeffRef(k,k);	327	A21 /= mat.coeffRef(k,k);
		328
		329	if(sign)
		330	{
		331	// LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
		332	int newSign = math::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
		333	if(k == 0)
		334	*sign = newSign;
		335	else if(*sign != newSign)
		336	*sign = 0;
		337	}
337	}	338	}
338		339
339	return true;	340	return true;
340	}	341	}
341		342
342	// Reference for the algorithm: Davis and Hager, "Multiple Rank	343	// Reference for the algorithm: Davis and Hager, "Multiple Rank
343	// Modifications of a Sparse Cholesky Factorization" (Algorithm 1)	344	// Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
344	// Trivial rearrangements of their computations (Timothy E. Holy)	345	// Trivial rearrangements of their computations (Timothy E. Holy)
Lines 362-391 template<> struct ldlt_inplace<Lower> Link Here
362	// Apply the update	363	// Apply the update
363	for (Index j = 0; j < size; j++)	364	for (Index j = 0; j < size; j++)
364	{	365	{
365	// Check for termination due to an original decomposition of low-rank	366	// Check for termination due to an original decomposition of low-rank
366	if (!(isfinite)(alpha))	367	if (!(isfinite)(alpha))
367	break;	368	break;
368		369
369	// Update the diagonal terms	370	// Update the diagonal terms
370	RealScalar dj = real(mat.coeff(j,j));	371	RealScalar dj = math::real(mat.coeff(j,j));
371	Scalar wj = w.coeff(j);	372	Scalar wj = w.coeff(j);
372	RealScalar swj2 = sigma*abs2(wj);	373	RealScalar swj2 = sigma*abs2(wj);
373	RealScalar gamma = dj*alpha + swj2;	374	RealScalar gamma = dj*alpha + swj2;
374		375
375	mat.coeffRef(j,j) += swj2/alpha;	376	mat.coeffRef(j,j) += swj2/alpha;
376	alpha += swj2/dj;	377	alpha += swj2/dj;
377		378
378		379
379	// Update the terms of L	380	// Update the terms of L
380	Index rs = size-j-1;	381	Index rs = size-j-1;
381	w.tail(rs) -= wj * mat.col(j).tail(rs);	382	w.tail(rs) -= wj * mat.col(j).tail(rs);
382	if(gamma != 0)	383	if(gamma != 0)
383	mat.col(j).tail(rs) += (sigmaconj(wj)/gamma)w.tail(rs);	384	mat.col(j).tail(rs) += (sigmamath::conj(wj)/gamma)w.tail(rs);
384	}	385	}
385	return true;	386	return true;
386	}	387	}
387		388
388	template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>	389	template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
389	static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)	390	static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
390	{	391	{
391	// Apply the permutation to the input w	392	// Apply the permutation to the input w

Lines 227-243 static typename MatrixType::Index llt_ra Link Here

(-)a/Eigen/src/Cholesky/LLT.h (-3 / +3 lines)
227	}	227	}
228	}	228	}
229	else	229	else
230	{	230	{
231	temp = vec;	231	temp = vec;
232	RealScalar beta = 1;	232	RealScalar beta = 1;
233	for(Index j=0; j<n; ++j)	233	for(Index j=0; j<n; ++j)
234	{	234	{
235	RealScalar Ljj = real(mat.coeff(j,j));	235	RealScalar Ljj = math::real(mat.coeff(j,j));
236	RealScalar dj = abs2(Ljj);	236	RealScalar dj = abs2(Ljj);
237	Scalar wj = temp.coeff(j);	237	Scalar wj = temp.coeff(j);
238	RealScalar swj2 = sigma*abs2(wj);	238	RealScalar swj2 = sigma*abs2(wj);
239	RealScalar gamma = dj*beta + swj2;	239	RealScalar gamma = dj*beta + swj2;
240		240
241	RealScalar x = dj + swj2/beta;	241	RealScalar x = dj + swj2/beta;
242	if (x<=RealScalar(0))	242	if (x<=RealScalar(0))
243	return j;	243	return j;
Lines 246-262 static typename MatrixType::Index llt_ra Link Here
246	beta += swj2/dj;	246	beta += swj2/dj;
247		247
248	// Update the terms of L	248	// Update the terms of L
249	Index rs = n-j-1;	249	Index rs = n-j-1;
250	if(rs)	250	if(rs)
251	{	251	{
252	temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);	252	temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
253	if(gamma != 0)	253	if(gamma != 0)
254	mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigmaconj(wj)/gamma)temp.tail(rs);	254	mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigmamath::conj(wj)/gamma)temp.tail(rs);
255	}	255	}
256	}	256	}
257	}	257	}
258	return -1;	258	return -1;
259	}	259	}
260		260
261	template<typename Scalar> struct llt_inplace<Scalar, Lower>	261	template<typename Scalar> struct llt_inplace<Scalar, Lower>
262	{	262	{
Lines 272-288 template<typename Scalar> struct llt_inp Link Here
272	for(Index k = 0; k < size; ++k)	272	for(Index k = 0; k < size; ++k)
273	{	273	{
274	Index rs = size-k-1; // remaining size	274	Index rs = size-k-1; // remaining size
275		275
276	Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);	276	Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
277	Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);	277	Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
278	Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);	278	Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
279		279
280	RealScalar x = real(mat.coeff(k,k));	280	RealScalar x = math::real(mat.coeff(k,k));
281	if (k>0) x -= A10.squaredNorm();	281	if (k>0) x -= A10.squaredNorm();
282	if (x<=RealScalar(0))	282	if (x<=RealScalar(0))
283	return k;	283	return k;
284	mat.coeffRef(k,k) = x = sqrt(x);	284	mat.coeffRef(k,k) = x = sqrt(x);
285	if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();	285	if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
286	if (rs>0) A21 *= RealScalar(1)/x;	286	if (rs>0) A21 *= RealScalar(1)/x;
287	}	287	}
288	return -1;	288	return -1;

Lines 263-282 template<typename MatrixType> void chole Link Here

(-)a/test/cholesky.cpp (-4 / +12 lines)
263		263
264	// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.	264	// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
265	// This test checks that LDLT reports correctly that matrix is indefinite.	265	// This test checks that LDLT reports correctly that matrix is indefinite.
266	// See http://forum.kde.org/viewtopic.php?f=74&t=106942	266	// See http://forum.kde.org/viewtopic.php?f=74&t=106942
267	template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)	267	template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
268	{	268	{
269	eigen_assert(m.rows() == 2 && m.cols() == 2);	269	eigen_assert(m.rows() == 2 && m.cols() == 2);
270	MatrixType mat;	270	MatrixType mat;
271	mat << 1, 0, 0, -1;	271	{
272	LDLT<MatrixType> ldlt(mat);	272	mat << 1, 0, 0, -1;
273	VERIFY(!ldlt.isNegative());	273	LDLT<MatrixType> ldlt(mat);
274	VERIFY(!ldlt.isPositive());	274	VERIFY(!ldlt.isNegative());
		275	VERIFY(!ldlt.isPositive());
		276	}
		277	{
		278	mat << 1, 2, 2, 1;
		279	LDLT<MatrixType> ldlt(mat);
		280	VERIFY(!ldlt.isNegative());
		281	VERIFY(!ldlt.isPositive());
		282	}
275	}	283	}
276		284
277	template<typename MatrixType> void cholesky_verify_assert()	285	template<typename MatrixType> void cholesky_verify_assert()
278	{	286	{
279	MatrixType tmp;	287	MatrixType tmp;
280		288
281	LLT<MatrixType> llt;	289	LLT<MatrixType> llt;
282	VERIFY_RAISES_ASSERT(llt.matrixL())	290	VERIFY_RAISES_ASSERT(llt.matrixL())