Cholesky modifications for A'A

chol_ata_append_col.m

function [l1, l2, flag]=chol_ata_append_col(L, A, c)
%CHOL_ATA_APPEND_COL updates the cholesky factorization of matrix A'A when a
%column is appended at the end of A, i.e., 
%
%knowning that the cholesky factorization of A'A is
%
% A'A = LL'
%
%we can efficiently compute the Cholesky factorization of the matrix
%
% [A c]'[A c]
%
%provided that it exists. This, will then have the form
%
%  [L   0
%  l1' l2]
%
%The total number of flops needed to perform this update is n^2+3n.
%
%Syntax:
%[l1, l2, flag]=CHOL_ATA_APPEND_COL(L, A, c)
%
%Input arguments:
%c  The column which is appended at the end of matrix A
%A  The original matrix A
%L  The Cholesky factor of A'A
%
%Output arguments:
%l1, l2    vectors which are used to update the Cholesky factorization of 
%          A'A when a column is appended in A
%flag      this flag is set to 0 if the computation has succeded and to 
%          -1 if the new matrix [A c]'[A c] is not positive definite.
%
%See also:
% chol_ata_insert_col, chol_ata_remove_col

% Pantelis Sopasakis

narginchk(3,3);
l1 = L\(A'*c);
l2 = c'*c - l1'*l1;
if l2>0,
    flag = 0;
    l2 = sqrt(l2);
else
    flag = -1;
    l2 = [];
end

chol_ata_insert_col.m

function [l1, l2, P_tr, flag]=chol_ata_insert_col(L, A, c, idx)
%CHOL_ATA_INSERT_COL is similar to chol_ata_append_col, but the difference
%is that here the additional column is inserted in A. Let A be a matrix
%in the following form
%
% A = [a(1) ... a(idx-1) a(idx) ... a(n)]
%
%where a(i) are column-vectors. We insert a column-vector c to get the
%matrix
%
% A_ = [a(1) ... a(idx-1) c a(idx) ... a(n)]
%
%This function computes the permuted Cholesky factorisation of A_'*A_ given
%the Cholesky factorisation of A'*A, that is, it constructs a lower
%triangular matrix L_ and a permutation matrix P_ so that
%
% A_'A_ = P_ L_ L_' P_'
%
%The matrix L_ has the following form
%
% L_ = [ L    0
%        l1'  l2]
%
%The (transpose of the) permutation matrix P_' is returned as a vector.
%Notice that P_' is such that A_*P_' = [A c]. Using the returned vector of
%integer indices P_, the following holds true:
%
% Ac'*Ac == T(P_tr, P_tr))
% 
%Syntax:
% [l1, l2, P_tr, flag]=CHOL_ATA_INSERT_COL(L, A, c, idx)
%
%Input arguments:
% L         Lower triangular matrix. The Cholesky factorisation of A'*A.
% A         Original matrix A
% c         Column to add in A
% idx       Index where c will be added
%
%Output arguments:
% l1, l2    Data used to update the Cholesky factorisation of A'A
% P_tr      Transpose of the permutation matrix P as described above
% flag      this flag is set to 0 if the computation has succeded and to 
%          -1 if the new matrix [A c]'[A c] is not positive definite.
%
%See also:
% chol_ata_append_col, chol_ata_remove_col



% Pantelis Sopasakis

narginchk(4,4);
n = size(A, 2);
P_tr = [1:idx-1 n+1 idx:n];
[l1, l2, flag] = chol_ata_append_col(L,A,c);

chol_ata_remove_col.m

function [L11bar, L31bar, L33bar] = chol_ata_remove_col(L, idx)
%CHOL_ATA_REMOVE_COL updated the Cholesky factorisation of A'A when a
%column is removed from matrix A. We assume that A has the form
%
% A = [a(1) a(2) ... a(idx-1) a(idx) a(idx+1) ... a(n)]
% 
%where a(i) are column-vectors and a(idx) is removed from A to form A_.
%Here is an example in MATLAB code:
%
% > A = rand(140, 12);
% > idx = 5;
% > A_ = [A(:,1:idx-1) A(:,idx+1:end)];
%
%The resulting matrix A_ leads to the following factorisation of A_'A_
%
% A_'A_ = [ L11bar  0
%           L31bar  L33bar ];
%
%This function computes matrices L11bar, L31bar and L33bar given the
%Cholesky factorisation of the original matrix.
%
%The total number of flops to perform this update is 2(n-idx)^2+4(n-idx).
%
%Syntax
%[L11bar, L31bar, L33bar] = CHOL_ATA_REMOVE_COL(L, idx)
%
%Input arguments:
% L     Cholesky factorisation of A'*A (lower triangular)
% idx   Integer index or indices of the columns in A to be removed
%
%
%Example:
% 
%  A = randn(150,12);
%  L = chol(A'*A,'lower');
%  cols_to_remove = [3 1 5 10];
%  [L11bar, L31bar, L33bar] = chol_ata_remove_col(L, cols_to_remove);
%  L_updated= [L11bar      zeros(6, 2);
%              L31bar      L33bar];
%
%
%See also:
% chol_ata_append_col, chol_ata_insert_col

% Pantelis Sopasakis

narginchk(2,2);

idx = sort(idx);

if (length(idx) > 1),
    L_updated = L;
    idx_rem = 0;
    for i=idx, 
        [L11bar, L31bar, L33bar] = chol_3587(L_updated, i-idx_rem);
        L_updated= [L11bar      zeros(size(L11bar,1), size(L33bar, 2));
                    L31bar      L33bar];
        idx_rem = idx_rem + 1;
    end
else
    [L11bar, L31bar, L33bar] = chol_3587(L, idx);
end

function [L11bar, L31bar, L33bar] = chol_3587(L, idx)
    L11bar=L(1:idx-1, 1:idx-1);
    l32 = L(idx+1:end, idx);
    L31bar=L(idx+1:end, 1:idx-1);
    L33bar=cholupdate(L(idx+1:end, idx+1:end)', l32, '+')';

test_chol_ata_append.m

clc;
clear;


%% UPDATE WHEN A NEW COLUMN IS APPENDED AT THE END

A=randn(150,7);               % random rectangular matrix
L = chol(A'*A,'lower');     % cholesky of original matrix

c=rand(size(A,1), 1);       % column to add

Ac = [A c];                 % new matrix

% compute new Cholesky (provide only L, A and c)
[l1, l2, flag] = chol_ata_append_col(L, A, c);
assert(flag==0);
assert(~isempty(l2));
assert(~isempty(l1));

% Construct the updated Cholesky
L_update = [L zeros(size(L,1),1); l1' l2];

% test
assert(norm(Ac'*Ac - L_update*L_update')<1e-10);

%% UPDATE WHEN A ROW IS APPENDED
A=randn(20,7);
L = chol(A'*A);
w = rand(size(A,2),1);
L_row = cholupdate(L, w,'+')';
assert(norm([A; w']'*[A; w']-L_row*L_row')<1e-10)

%%
clc
N=5000;
A=rand(N,100);
L=chol(A'*A,'lower');

c=rand(N,1);
[l1, l2,flag] = chol_ata_append_col(L,A,c);
L_update = [L zeros(size(L,1),1); l1' l2];
assert(flag==0);

d=c+0.00001*randn(N,1);
[l1, l2, flag] = chol_ata_append_col(L_update,[A c],d);;
assert(flag==0);

test_chol_ata_insert.m

clc;
clear
A = randn(150,7);
L = chol(A'*A,'lower');

c=rand(size(A,1), 1); 
idx = 4;


Ac = [A(:,1:idx-1) c A(:, idx:end) ];

[l1, l2, P_tr, flag]=chol_ata_insert_col(L, A, c, idx);
[~, P_] = sort(P_tr);
assert(flag==0);

assert( norm([A c]-Ac(:,P_)) < 1e-14  ); % make sure the permutation is correct

I=eye(8);
Pmat = I(:, P_tr);


L_updated = [L zeros(size(L,1),1); l1' l2];

T=L_updated*L_updated';
assert(  norm(T-[A c]'*[A c]) < 1e-10  );
assert(  norm(Ac'*Ac - Pmat'*T*Pmat) < 1e-10  );
assert(  norm(Ac'*Ac - T(P_tr, P_tr)) < 1e-10  );

test_chol_ata_remove.m

% ! rm -r *~ ./tests/*~
clc;
clear;

for r=1:100,
    A = randn(150,7);                 % random rectangular matrix
    L = chol(A'*A,'lower');           % cholesky of original matrix


    for idx = 1:7,
        [L11bar, L31bar, L33bar] = chol_ata_remove_col(L, idx);

        L_updated= [L11bar      zeros(idx-1,size(A,2)-idx);
                    L31bar      L33bar];

        A_removed_col = [A(:,1:idx-1) A(:,idx+1:end)];
        L_correct = chol(A_removed_col'*A_removed_col,'lower');

        assert(norm(L_updated-L_correct)<1e-10);
    end
end

%% Remove many columns
clear
for r=1:100,
    A = randn(150,12);
    L = chol(A'*A,'lower');
    cols_to_remove = [1 3 5 10];
    cols_to_remove = sort(cols_to_remove);
    sd = setdiff((1:12), cols_to_remove);
    Anew = A(:, sd);
    L_updated = L;
    idx_rem = 0;
    for idx=cols_to_remove,   
        [L11bar, L31bar, L33bar] = chol_ata_remove_col(L_updated, idx-idx_rem);
        L_updated= [L11bar      zeros(size(L11bar,1), size(L33bar, 2));
                    L31bar      L33bar];               
        idx_rem = idx_rem+1;
    end

    assert(    norm(Anew'*Anew - L_updated*L_updated')   < 1e-10 );
end

%% Remove many one command
clear

for r=1:100,
    A = randn(150,12);
    L = chol(A'*A,'lower');
    cols_to_remove = [3 1 5 10];
    [L11bar, L31bar, L33bar] = chol_ata_remove_col(L, cols_to_remove);
    L_updated= [L11bar      zeros(6, 2);
        L31bar      L33bar];


    sd = setdiff((1:12), cols_to_remove);
    Anew = A(:, sd);   
    assert(    norm(Anew'*Anew - L_updated*L_updated')   < 1e-10 );
end