cholrank1 update with LDL decomposition

97 Views Asked by At

I have cholrank1 update procedure (wikipedia) for the symmetric positive definite (SPD) matrix.

function [L] = cholupdate(L,x)
    p = length(x);
    for k=1:p
        r = sqrt(L(k,k)^2 + x(k)^2);
        c = r / L(k, k);
        s = x(k) / L(k, k);
        L(k, k) = r;
        L(k+1:p,k) = (L(k+1:p,k) + s*x(k+1:p)) / c;
        x(k+1:p) = c*x(k+1:p) - s*L(k+1:p,k);
    end
end

It works with LL decomposition. I try to fix procedure to work with LDL decomposition (ie without calling sqrt) like this:

function [L] = cholupdate_ldl(L,x)
    p = length(x);
    for k=1:p
        r = L(k,k) + x(k)^2;
        c = r / L(k, k);
        s = x(k) / L(k, k);
        L(k, k) = r;
        L(k+1:p,k) = (L(k+1:p,k) + s*x(k+1:p)) / c;
        x(k+1:p) = sqrt(c)*(x(k+1:p) - x(k)*L(k+1:p,k));
    end
end

It works fine but I was forced to use sqrt. How can I update LDL decomposition without using sqrt at all?