## How to dynamically construct a matrix of functions ?

Asked by Thomas Barbier

### Thomas Barbier (view profile)

on 20 May 2019
Latest activity Commented on by Fangjun Jiang

### Fangjun Jiang (view profile)

on 24 May 2019
Dear all,
I do have two functions and that both return a matrix (let's say 2x2 for the sake of the argument). I would like to construct a biggerA matrix function which takes some arguments such that:
I cant "hard-code" the handles as the size of A may not always be the same: if and are available thenA would be:
Is this possible?
Thanks for your help.
Regards

R2018b

### Fangjun Jiang (view profile)

Answer by Fangjun Jiang

### Fangjun Jiang (view profile)

on 20 May 2019
Edited by Fangjun Jiang

### Fangjun Jiang (view profile)

on 20 May 2019

I think you can make the matrix using
M1=vertcat(blkdiag(JcX1,JcX2), blkdiag(JcX3,JcX4), ...)
M2=blkdiag(vertcat(JfX1,JfX2), vertcat(JfX1,JfX2), ...)
A=horzcat(M1,M2)

Thomas Barbier

### Thomas Barbier (view profile)

on 23 May 2019
Yes that's it. Sorry the problem is already not easy to grasp and my presentation of doesn't help.
At the time of writing the code, the size of x is unknown. The size of x depends of the problem at the time of the run. On some runs x may be 100, on some others x may be 300.
When the code is run, x is known : it's a first estimate. Thus the shape of A is known and does not change, but it's numerical value does change along the optimization (the loop in "nonLinearRegression").
Hence my quest to construct at the beginning of the run (so I know the size of x) a "symbolical" matrix @(x) A(x) that I can then call during the optimisation (with the numerical values of x_num). That would avoid reconstructing a numerical matrix at each loop of optimization and just doing A_num = A(x_num).
Thank you
Fangjun Jiang

### Fangjun Jiang (view profile)

on 24 May 2019
This should do it.
x=1:6;
Jc=@(u) u-0.1+ones(2);
Jf=@(u) u-0.9+ones(2);
A=ZeroA(x);
for k=1:numel(x)
A=FillA(A,k,x(k),Jc,Jf);
end
open A;
function A=ZeroA(x)
n=numel(x);
A=zeros(2*n, 4+n);
end
function A=FillA(A,n,Xn,Jc,Jf)
IndX=(-1:0)+2*n;
IndY=IndX-4*ceil(n/2-1);
A(IndX,IndY)=Jc(Xn);
IndY=(3:4)+2*ceil(n/2);
A(IndX,IndY)=Jf(Xn);
end
Fangjun Jiang

### Fangjun Jiang (view profile)

on 24 May 2019
Function ZeroA() creates the matrix A with the proper size based on X=[X1, X2, ..., Xn].
Function FillA() fills the proper elments of A based on one element of X (for example X2). It does provide efficiency if X is 1x300 and chenges only one element at a time.
Any way, it seems much harder to understand a question than to anwer one.

### Guillaume (view profile)

on 23 May 2019

If I understood correctly:
function A = makeMatrix(Jc, Jf, x)
%Jc: function handle to a function that takes scalar input xi and returns a fixed size matrix
%Jf: function handle to a function that takes scalar input xi and returns a fixed size matrix
%x: vector of xi, must have an EVEN number of elements
assert(mod(numel(x), 2) == 0, 'Precondition broken');
Jcx = arrayfun(Jc, x, 'UniformOutput', false); %apply Jc to each x, store result in cell array
Jfx = arrayfun(Jf, x, 'UniformOutput', false); %apply Jf to each x, store result in cell array
Jcz = repmat({zeros(size(Jcx{1}))}, 1, numel(x)/2); %matrices of zeros to insert in first two columns
Jcmat = cell2mat(reshape([Jcx{1:2:end}, Jcz; Jcz, Jcx{2:2:end}], [], 2)); %build first two columns
Jfx = cellfun(@(odd, even) [odd;even], Jfx(1:2:end), Jfx(2:2:end), 'UniformOutput', false); %concatenate odd and even consecutive indices
Jfmat = blkdiag(Jfx{:}); %blkdiag the odd-even matrices
A = [Jcmat, Jfmat];
end

Show 1 older comment
Guillaume

### Guillaume (view profile)

on 24 May 2019
The only way you could do this is with the symbolic toolbox. Unfortunately, I don't know anything about the symbolic toolbox.
I don't see how the accepted answer does what you want. It also reconstructs A at each call.
Thomas Barbier

### Thomas Barbier (view profile)

on 24 May 2019
Yes indeed, I came with the same conclusion. So I'll reconstruct each time...
Guillaume

### Guillaume (view profile)

on 24 May 2019
Without the symbolic toolbox, what you could do to speed up the reconstruction is store the location of the Jcx and Jfx in the final matrix beforehand, so the only thing you have to do in the loop is calculate the actual values and insert them at the precalculated locations:
function [locJc, locJf] = PrepareA(xlength, matsize)
%xlength: number of elements in x. Must be even
%matsize: two element vectors indicated the size of matrices returned by Jcx and Jfx
%locJx: 2 columns matrices indicating the insertion point of Jxx
locJc = [(0:xlength-1)*matsize(1) + 1; repmat([1, matsize(2)+1], 1, xlength/2)].';
locJf = [(0:xlength-1)*matsize(1) + 1; repelem(0:xlength/2-1, 2)*matsize(2) + 2*matsize(2) + 1].';
end
When it's time to reconstruct A:
function A = constructA(Jc, Jf, x, locJc, locJf, matsize)
%Jc, Jf: function handles
%x: array of x values
A = zeros(matsize(1) * numel(x), matsize(2) * (numel(x)/2 + 2));
for xidx = 1:numel(x)
A(locJc(xidx, 1) + (0:matsize(1)-1), locJc(xidx, 2) + (0:matsize(2)-1)) = Jc(x(xidx));
A(locJf(xidx, 1) + (0:matsize(1)-1), locJf(xidx, 2) + (0:matsize(2)-1)) = Jf(x(xidx));
end
end
To speed up even more, you could make locJc and locJf cell arrays of indices so as to avoid any index calculation in constructA.