Asked by Thomas Barbier
on 20 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

Answer by Fangjun Jiang
on 20 May 2019

Edited by Fangjun Jiang
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
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
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
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.

Sign in to comment.

Answer by Guillaume
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

Guillaume
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
on 24 May 2019

Yes indeed, I came with the same conclusion. So I'll reconstruct each time...

Guillaume
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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.