Concatenate along arbitrary dimension

[Wannes Van Loock]

At the moment there is only support for horizontal or vertical concatenation of sdpvars. Concatenating along another dimension is not possible. Is this feature planned in a future release?

[Johan Löfberg]

That is a fair request, so I'll add it to the list.

Speaking of which, I guess I should fix this bug also

>> x = sdpvar(3,3,3);
>> y = sdpvar(3,3,3);
>> [x y]
Multi-dimensional SDPVAR object 6x3x3
>> [x;y]
Multi-dimensional SDPVAR object 6x3x3

[Johan Löfberg]

Place these in extras/@ndsdpvar

[Wannes Van Loock]

Wow, thanks for the lightning fast support!

For ndsdpvar it's working, but cat should still be overloaded for regular sdpvar objects. I tried converting sdpvar to ndsdpvar but then cat is complaining about the sizes...

[Johan Löfberg]

You mean like this

>> x = randn(3,3,3);
>> y = randn(3,3);
>> z=cat(3,x,y);
>> w=cat(3,y,y)

[Wannes Van Loock]

Exactly.

[Johan Löfberg]

As a temporary fix, replace the dimension check in cat.m with

m = max([length(dimx) length(dimy) along]);
dimx = [dimx ones(1,m-length(dimx))];
dimy = [dimy ones(1,m-length(dimy))];

and the cast the sdpvar as an ndsdpvar when you want to do concatenation to higher dimensions

>> x = sdpvar(2,2);
>> y = sdpvar(2,2);
>> cat(3,x,y)
>> cat(3,ndsdpvar(x),ndsdpvar(y))
Multi-dimensional SDPVAR object 2x2x2

[Johan Löfberg]

change sdpvar/cat

function y = cat(varargin)
%CAT (overloaded)

switch varargin{1}
    case 1
        y = vertcat(varargin{2:end});
    case 2
        y = horzcat(varargin{2:end});
    otherwise
        for i = 2:nargin
            if isa(varargin{i},'sdpvar');
                varargin{i} = ndsdpvar(varargin{i});
            end
        end
        y = cat(varargin{:});
end

and you don't have to cast manually

[Wannes Van Loock]

Thanks for the quick fix!

Related to this issue: ndims should also be overloaded. I simply added to the ndsdpvar class

function n = ndims(X)
    n = length(size(X))

[Wannes Van Loock]

Performancewise this implementation doesn't seem optimal.

In my application, 97% (140s for a very small problem) is spent in ndsdpvar/cat, where most time is lost in the lines.

z = sdpvar(length(zindex),1);
z(find(locx))=sdpvar(x(:));
z(find(locy))=sdpvar(y(:));

As an alternative I could try to make my code less dependent on cat.

[Johan Löfberg]

I was secretely hoping you had a very small problem. No, it is very inefficient so this is just a temporary fix. What dimensions are you targeting?

[Johan Löfberg]

Easily fixed. Replace those lines by

A = sparse(find(locx),1:prod(dimx),1,prod(zdim),prod(dimx));
B = sparse(find(locy),1:prod(dimy),1,prod(zdim),prod(dimy));
z = A*sdpvar(x(:)) + B*sdpvar(y(:));

[Wannes Van Loock]

I'm working with multivariate (nd) functions, where the function can be represented as a tensorproduct of some coefficients, the optimization variables, and basis functions (with dimensions k_1, k_2, ..., k_nd). The functions can be scalar, vector or matrix valued. These coefficients are stored in an (k_1 x k_2 x ... x k_nd)-cell, in which each element is an sdpvar (of dimension (m x m)). As you can see, these things quickly grow quite big. To do function manipulations, I require a tensor representation of the coefficients, for which I'm using a slightly modified version of cell2mat and this function relies heavily on cat.

Anyways, your proposed fix already gives a huge performance improvement!

[Johan Löfberg]

Are you saying cat still takes a significant part of the overall computations, despite the sparse matrix multiplication hack?

[Wannes Van Loock]

The code runs about a factor 10 times faster, but most time is still being spent in the concatenation. However, I should mention that ndsdpvar.cat is still called 1295 times.

[Johan Löfberg]

Would like to see or know more about the scenario. Are you calling cat with many arguments, or just two arrays to concatenate?

[Wannes Van Loock]

I'm calling cat with many arguments. Below you can find an example for a matrix valued (3x3) 4-D multivariate functions, with functions basis of size (5x5x5x5).
As said before, most of this code is a stripped down version of cell2mat. For double arrays this is fast but perhaps the functionality can be implemented more efficiently for YALMIP variables.

c = arrayfun(@(~)sdpvar(3,3), zeros(5,5,5,5), 'uni', false);  % Create coefficient matrix
elements = numel(c);

if elements == 1
    m = c{1};
    return
end

csize = size(c);
% Construct the matrix by concatenating each dimension of the cell array into
%   a temporary cell array, CT
% The exterior loop iterates one time less than the number of dimensions,
%   and the final dimension (dimension 1) concatenation occurs after the loops

% Loop through the cell array dimensions in reverse order to perform the
%   sequential concatenations
for cdim=(length(csize)-1):-1:1
    % Pre-calculated outside the next loop for efficiency
    ct = cell([csize(1:cdim) 1]);
    cts = size(ct);
    ctsl = length(cts);
    mref = {};

    % Concatenate the dimension, (CDIM+1), at each element in the temporary cell
    %   array, CT
    for mind=1:prod(cts)
        [mref{1:ctsl}] = ind2sub(cts,mind);
        % Treat a size [N 1] array as size [N], since this is how the indices
        %   are found to calculate CT
        if ctsl==2 && cts(2)==1
            mref = {mref{1}};
        end
        % Perform the concatenation along the (CDIM+1) dimension
        ct{mref{:}} = cat(cdim+1,c{mref{:},:});
    end
    % Replace M with the new temporarily concatenated cell array, CT
    c = ct;
end
% Finally, concatenate the final rows of cells into a matrix
m = cat(1,c{:});

[Johan Löfberg]

The attached file is a bit faster (but will require some more coding before being fully functional for more complex scenarios. Works for your cases though)

[Johan Löfberg]

BTW, I think you can cut off some time with

c = sdpvar(3*ones(1,625),3*ones(1,625));
c = reshape(c,[5 5 5 5]);