Permalink
Browse files

slightly improved handling of singular JtJ

Previously I was simply adding a constnat lambda*I to JtJ if it was ever
singular. I'm now still doing this, but at least I increase lambda every time
it's singular. I never reduce it. And this is still a somewhat silly method.
  • Loading branch information...
dkogan committed Mar 25, 2012
1 parent 36de44c commit c36211688eb66eb5587eaa66227f89a2989e404a
Showing with 32 additions and 23 deletions.
  1. +27 −19 dogleg.c
  2. +5 −4 dogleg.h
View
@@ -44,7 +44,7 @@ static double UPDATE_THRESHOLD = 1e-8;
static double TRUSTREGION_THRESHOLD = 1e-8;
// if I ever see a singular JtJ, I factor JtJ + LAMBDA*I from that point on
-#define LAMBDA 1e-20
+#define LAMBDA_INITIAL 1e-10
// these parameters likely should be messed with
void dogleg_setDebug(int debug)
@@ -296,23 +296,31 @@ static void computeGaussNewtonUpdate(dogleg_operatingPoint_t* point, dogleg_solv
}
// try to factorize the matrix directly. If it's singular, add a small
- // constant to the diagonal from this point on
- if(!ctx->wasPositiveSemidefinite)
+ // constant to the diagonal. This constant gets larger if we keep being
+ // singular
+ while(1)
{
- ASSERT( cholmod_factorize(point->Jt, ctx->factorization, &ctx->common) );
- if(ctx->factorization->minor != ctx->factorization->n)
+ if( ctx->lambda == 0.0 )
+ ASSERT( cholmod_factorize(point->Jt, ctx->factorization, &ctx->common) );
+ else
{
- ctx->wasPositiveSemidefinite = 1;
- if( DOGLEG_DEBUG )
- fprintf(stderr, "singular JtJ. Adding %g I from now on\n", LAMBDA);
+ double beta[] = { ctx->lambda, 0 };
+ ASSERT( cholmod_factorize_p(point->Jt, beta, NULL, 0,
+ ctx->factorization, &ctx->common) );
}
+
+ if(ctx->factorization->minor == ctx->factorization->n)
+ break;
+
+ // singular JtJ. Raise lambda and go again
+ if( ctx->lambda == 0.0) ctx->lambda = LAMBDA_INITIAL;
+ else ctx->lambda *= 10.0;
+
+ if( DOGLEG_DEBUG )
+ fprintf(stderr, "singular JtJ. Have rank/full rank: %zd/%zd. Adding %g I from now on\n",
+ ctx->factorization->minor, ctx->factorization->n, ctx->lambda);
}
- if(ctx->wasPositiveSemidefinite)
- {
- double beta[] = {LAMBDA, 0};
- ASSERT( cholmod_factorize_p(point->Jt, beta, NULL, 0,
- ctx->factorization, &ctx->common) );
- }
+
// solve JtJ*updateGN = Jt*x. Gauss-Newton step is then -updateGN
cholmod_dense Jt_x_dense = {.nrow = point->Jt->nrow,
@@ -705,11 +713,11 @@ double dogleg_optimize(double* p, unsigned int Nstate,
void* cookie,
dogleg_solverContext_t** returnContext)
{
- dogleg_solverContext_t* ctx = malloc(sizeof(dogleg_solverContext_t));
- ctx->f = f;
- ctx->cookie = cookie;
- ctx->factorization = NULL;
- ctx->wasPositiveSemidefinite = 0;
+ dogleg_solverContext_t* ctx = malloc(sizeof(dogleg_solverContext_t));
+ ctx->f = f;
+ ctx->cookie = cookie;
+ ctx->factorization = NULL;
+ ctx->lambda = 0.0;
if( returnContext != NULL )
*returnContext = ctx;
View
@@ -50,10 +50,11 @@ typedef struct
// the factorization of the most recent JtJ
cholmod_factor* factorization;
- // Have I ever seen a singular JtJ? If so, I add a small constant to the
- // diagonal from that point on. This is a simple and fast way to deal with
- // singularities. This is suboptimal but works for me for now.
- int wasPositiveSemidefinite;
+ // Have I ever seen a singular JtJ? If so, I add this constant to the diagonal
+ // from that point on. This is a simple and fast way to deal with
+ // singularities. This constant starts at 0, and is increased every time a
+ // singular JtJ is encountered. This is suboptimal but works for me for now.
+ double lambda;
} dogleg_solverContext_t;

0 comments on commit c362116

Please sign in to comment.