How to dynamically construct a matrix of functions ?

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 1

[Guillaume]

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

[Thomas Barbier]

Thank you for your answer Guillaume.

That would work except that it "reconstructs" the numericlal matrix at each call. Infortunately that's not what I aim for.

When I run the code, the size of is then known. From that, I'd like to construct a "symbolical" version of (with the pattern shown in the first question) so that I can evaluate a numerical version of whenever needed just by doing:

A_num = A(x_num);

[Guillaume]

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]

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

[Guillaume]

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.

Answer 2

[Fangjun Jiang]

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]

Hi,

Thank you for your answer. This would work if I knew in advance the final size of my matrix A. However it's size depends of the size of which I know only when I want to have the numerical expression of the matrix.

To explain a bit more, this is part of an optimization scheme, where the matrix is the Jacobian of the state vector . The optimization has a state vector with sometimes and sometimes with (depending on measurements).

The optimization is done thanks to a non-linear regression algorithm (such as Gauss-Newton or Levenberg-Marquardt). The the state vector is defined at the beginning of the problem, thus the size of the Matrix won't change during all the iterations of the regression algorithm. I could build "by hand" But the algorithm can iterate up to 100 times, and the numerical values change at each iteration. So I thought that defining a function with the Jacobian shape would be more efficient than assembling numerical blocks at each iteration.

That could be sumed up with:

function solution = optimize(x, f)
    % x the state-vector that needs to be optimized
    % f the function on which evolves x
    len = length(x);
    
    % Construct the matrix A which is dependant of the size of x
    
    solution = nonLinearRegression(x, f, A);
end

function xOptimized = nonLinearRegression(x, F, JacobianFunction)
    residualMin = 10e-6;
    
    do {
      %Computation of the numerical Jacobian at the operating point
      Jnum = JacobianFunction(x);
      %Computation of the numerical function at the operating point
      Fnum = F(x);
      
      Delta = - inv(Jnum' * Jnum) * Jnum' * Fnum;
      x = x + Delta;
    } while (Delta > resildualMin)

end

[Fangjun Jiang]

I don't follow.

The three lines of code didn't pre-specify the size of A. The size of A is the result of blkdiag(), vercat() and horzcat().

[Thomas Barbier]

Before executing the code, I don't know the size of , which prevents me using vertcat or blkdiag. I want to have a function at the end A(x).

I can't make:

len = length(x);

A = @(x) ...
for i = 1:n
    A(i+offset1, i+offset2) = Jf(x(i))
end

[Fangjun Jiang]

So, can we assume

• x is variable size, its size could be 1x6, 1x8, 1x10, etc. ?
• Jc(x1), Jf(x1), Jc(x2), Jf(x2), etc. are all scalars ?

[Guillaume]

I understood that Jc(x1), etc. are matrices "I do have two functions Jc(x) and Jf(x) that both return a matrix."

[Thomas Barbier]

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]

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]

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.