Error using ichol: Encountered nonpositive pivot

[economoser]

First of all, thank you so much for the fantastic work and sharing this with the community!

I’ve noticed that sometimes in one of the iterations the infamous "Error using ichol": "Encountered nonpositive pivot" can occur when computing the incomplete Cholesky factorization of the Gramian matrix as a preconditioner for matrix inversion in the estimation. For example, I believe this can happen in the call to the ichol function that occurs six lines after the comment "%%Deresidualized [...]" in version 2.1 of leave_out_COMPLETE.m. I believe -- TBC! -- I have also seen it occur later on in the code during another call to the same ichol function.

When this error occurs, then currently the whole run is aborted. I wonder if there’s a way to -try- and -catch- that error or otherwise find some workaround that lets you still use the information obtained during the previous code evaluations? For example, when there are many calls to ichol during the random projection iterations, could a failed iteration be replaced with another one? This would be a marvelous feature as sometimes the routine runs for a couple days before hitting an error.

NB: I am sure that my design matrix (X in y = X*beta + e) is not singular, i.e., I have dropped the necessary indicators and other covariates to avoid multicollinearity in controls.

[economoser]

Just confirming that the same error still occurs in code version 3.02:

-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
SECTION 2
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
{Error using ichol
Encountered nonpositive pivot.

Error in leave_out_KSS (line 450)
   Lchol=ichol(xx,struct('type','ict','droptol',1e-2,'diagcomp',.1));
}

[rsaggio87]

I thought a bit about this and I'm not sure what is the best actions. In this new version ichol is not required for computation of the leverages so this error should not impact the computation of the statistical leverages. I tried "catching" the error in the following way, let me know if you think it makes sense (this is in the Branch 3.1. waiting on your input to push this to master).

[economoser]

What you suggest is a very nice solution. If I may suggest one additional tweak: before resorting to brute force conjugate gradient method (PCG), one could loop through a number of diagcomp parameters to try more diagonally dominant matrices in the incomplete Cholesky decomposition, something along the lines of:

if no_controls == 0
    S=speye(J-1);
    S=[S;sparse(-zeros(1,J-1))];  %N+JxN+J-1 restriction matrix 
    X=[D,F*S,controls];
    xx=X'*X;
    xy=X'*y;
    diagcomp_linear_n = 1;
    diagcomp_linear_step = .1;
    diagcomp_linear_list = diagcomp_linear_step.*(1:diagcomp_linear_n);
    diagcomp_max = max(sum(abs(xx), 2)./diag(xx)) - 2; % value recommended by MATLAB to guarantee successful execution of -ichol()-, see https://www.mathworks.com/help/matlab/ref/ichol.html
    diagcomp_candidate_n = 20;
    diagcomp_candidate_base = 1.5;
    diagcomp_overshoot = 3;
    diagcomp_factor = 1./(diagcomp_candidate_base.^(diagcomp_candidate_n:-1:-diagcomp_overshoot));
    diagcomp_candidates = diagcomp_max.*diagcomp_factor;
    diagcomp_exponential_list = diagcomp_candidates(diagcomp_candidates > diagcomp_linear_step*diagcomp_linear_n);
    diagcomp_list = [diagcomp_linear_list diagcomp_exponential_list];
    Lchol = [];
    for diagcomp = diagcomp_list
        try
            Lchol=ichol(xx,struct('type','ict','droptol',droptol,'diagcomp',diagcomp));
            fprintf('\n')
            disp(['NOTE: function -ichol()- with diagcomp = ' num2str(diagcomp) ' succeeded!'])
            break % exit for loop after successful evaluation of -ichol()-
        catch
            disp(['USER WARNING: function -ichol()- with diagcomp = ' num2str(diagcomp) ' failed!'])
            if diagcomp == diagcomp_list(end)
                fprintf('\n')
                disp(['USER WARNING: function -ichol()- did not execute successfully for any value of diagcomp <= ' num2str(diagcomp)])
            end
        end
    end
    if size(Lchol,1) > 0 % if one of the -ichol()- evaluations succeeded, then use preconditioner
        b=pcg(xx,xy,1e-10,1000,Lchol,Lchol');
    else % else, run brute-force conjugate-gradient method
        b=pcg(xx,xy,1e-10,1000);
    end
    y=y-X(:,N+J:end)*b(N+J:end); %variance decomposition will be based on this residualized outcome.
end

Here, diagcomp_list is a list of diagcomp values to try before proceeding with the brute-force PCG. It turns out that, at least in my trials, these -ichol()- evaluations are pretty fast, so trying a bunch of values is pretty cheap computationally and could save a lot of time.

[rsaggio87]

amazing, thank you. Your version is much better. Gonna push this now.