Unexpected error with HODLR matrix mldivide

[sgsellan]

Hi!

My thanks for making this toolbox available. Extremely useful!

I was testing it out using a simple squared exponential kernel to build the matrices, but run into an unexpected error (sometimes!) when trying to solve a linear system with a matrix right-hand side.

% Set of points
rng(0);
num_points = 1000;
P = rand([num_points,2]);

% Build K3 the traditional way
disp("Building K3 the traditional way")
tic;
K3 = K3eval(1:num_points,1:num_points,P);
toc;

% Building K3 the hodlr way
disp("Building K3 the hodlr way")
tic;
hodlrK3 = hodlr('handle',@(i,j) K3eval(i,j,P),num_points,num_points);
toc;

% Test points
num_test = 5000;
Ptest = rand([num_test,2]);

disp("Build K2 the traditional way");
tic;
K2 = K2eval(1:num_test,1:num_points,Ptest,P);
toc;

disp("Build K2 the hodlr way");
tic;
hodlrK2 = hodlr('handle',@(i,j) K2eval(i,j,Ptest,P),num_test,num_points);
toc;

disp("Solve K3\K2 the traditional way")
tic;
sol = K3\K2';
toc;

disp("Solve K3\K2 the hodlr way")
tic;
hodlrsol = hodlrK3\(hodlrK2');
toc;

disp(['Error: ',num2str(max(max(abs(full(hodlrsol) - sol))))])

function v=K3eval(I,J,P)
    num_i = length(I);
    num_j = length(J);
    dim = size(P,2);
    rsq = zeros(num_i,num_j);
    for dd=1:dim
        rsq = rsq + repmat(P(I,dd),1,num_j).*repmat(P(J,dd)',num_i,1);
    end
    rsq = rsq + (repmat(reshape(I,[],1),1,num_j)==repmat(reshape(J,1,[]),num_i,1))*0.1;
    v = exp(-rsq);
    assert(size(v,1)==num_i);
    assert(size(v,2)==num_j);
end

function v=K2eval(I,J,Ptest,Ptrain)
    num_i = length(I);
    num_j = length(J);
    dim = size(Ptest,2);
    rsq = zeros(num_i,num_j);
    for dd=1:dim
        rsq = rsq + repmat(Ptest(I,dd),1,num_j).*repmat(Ptrain(J,dd)',num_i,1);
    end
    v = exp(-rsq);
    assert(size(v,1)==num_i);
    assert(size(v,2)==num_j);
end

I get the following error:

Error using  \ 
Matrix dimensions must agree.

Error in solve_upper_triangular (line 5)
        H.F = H1.F\H2.F;

Error in  \  (line 15)
        H = solve_upper_triangular(H1, H2);

Error in solve_lower_triangular (line 14)
            H.A11 = H1.A11\H2.A11;

Error in solve_lower_triangular (line 11)
            H.A11 = solve_lower_triangular(H1.A11, H2.A11);

Error in  \  (line 17)
        H = solve_lower_triangular(H1, H2);

Error in  \  (line 20)
        H = HU \ (HL\H2);

Error in hodlr_test
hodlrsol = hodlrK3\(hodlrK2');

Could you confirm that this is not the expected behaviour? If so, any idea how to fix this?

Thanks!

Silvia

[sgsellan]

Ah, digging a little deeper I see this may be related to the clustering being different in the two matrices I am multiplying (since one of them is rectangular). Perhaps mldivide should check for cluster consistency?

[robol]

Hi, and thanks for the report. My first guess would have been indeed incompatible clustering. Not all routines check for compatibility beforehand, it would probably be a good idea to add it.

In your specific case, once you build K2 you can use its cluster to build K3 as well, so that they are force to be the same. You can get row and columns clusters of an existing HODLR matrix with the cluster method, and use it by calling the constructor as follows:

H = hodlr('handle', ..., 'cluster', c);

If you specify a single cluster as above it is used for both rows and columns. If you want different clustering for rows and columns, you need to specify them as hodlr(..., 'cluster', rows, cols);.