Vectorizing equations for fsolve

Question

I have a system of 49 equations and 49 unknowns (\$x_{1},...,x_{49}\$). All equations have the following form.

\$x_{s}\sum_{r=1}^{R}\sum_{t=1}^{T}a_{s,t}b_{rs,t}=\sum_{r=1}^{R}\sum_{t=1}^{T}a_{r,t}b_{sr,t}x_{r}\$

where \$s = 1,...,49\$ with \$S = R = 50\$ and \$T = 30\$. The \$a\$'s, the \$b\$'s and \$x_{50}\$ are known parameters. Since all equations have the same form, I believe that vectorizing the problem is possible. I appreciate any suggestions of how to do this.

---

The code below solves a smaller system with \$S=R=T=3\$. The code runs. I describe how the parameters are arranged beneath the code.

    %% Define parameters
        T = 3; % Number of types (T)
        S = 3; % Number of senders (S)
        R = S; % Number of receivers (R), equals the number of senders
        A = [.1 .1 .1   % Columns are senders (S), rows are types (T)
             .3 .3 .6
             .6 .6 .3]; 
        B1 = [.9 .0 .1     
              .5 .1 .2     
              .8 .2 .1];   
        B2 = [.0 .7 .2     
              .2 .9 .0     
              .1 .8 .1];   
        B3 = [.1 .3 .7     
              .3 .0 .8     
              .1 .0 .8];   
        B = zeros(T,R,S); % Initialize
        B(:,1,:) = B1;    % Pages of B are receivers
        B(:,2,:) = B2;    % Rows of B are types
        B(:,3,:) = B3;    % Columns of B are senders

        % Starting values
        x0 = [36 26];
        % Call functions
        F = fsolve(@(x) myfun1(x,A,B),x0);

    %% Define Functions
        function F = myfun1(x,A,B) 
        F(1) =   x(1)*(sum(A(:,1).*B(:,2,1))  +  sum(A(:,1).*B(:,3,1))) ...
              -  x(2)*(sum(A(:,2).*B(:,1,2))) - (sum(A(:,3).*B(:,1,3)));
        F(2) =   x(2)*(sum(A(:,2).*B(:,1,2))  +  sum(A(:,2).*B(:,3,2))) ...
              -  x(1)*(sum(A(:,1).*B(:,2,1))) - (sum(A(:,3).*B(:,2,3)));    
        end

1. Is it possible to avoid the "sum"? Maybe by rearranging the parameters?
2. Is it possible to vectorize the function itself to avoid writing F(1), F(2), ... , F(49)? I know it is possible, but for this particular problem I have difficulties understanding how to do it.

A visualization of the B array: In the example, the parameters are arranged as shown below.

\$A:\$ \begin{pmatrix} a_{1,1} & a_{2,1} & a_{3,1}\\ a_{1,2} & a_{2,2} & a_{3,2}\\ a_{1,3} & a_{2,3} & a_{3,3} \end{pmatrix}

\$B1:\$ \begin{pmatrix} b_{11,1} & b_{12,1} & b_{13,1}\\ b_{21,1} & b_{22,1} & b_{23,1}\\ b_{31,1} & b_{32,1} & b_{33,1} \end{pmatrix}

\$B2:\$ \begin{pmatrix} b_{11,2} & b_{12,2} & b_{13,2}\\ b_{21,2} & b_{22,2} & b_{23,2}\\ b_{31,2} & b_{32,2} & b_{33,2} \end{pmatrix}

\$B3:\$ \begin{pmatrix} b_{11,3} & b_{12,3} & b_{13,3}\\ b_{21,3} & b_{22,3} & b_{23,3}\\ b_{31,3} & b_{32,3} & b_{33,3} \end{pmatrix}

Answer 1

Yes. You can simplify this code. I made two main changes.

First, I forget that \$x_{50}\$ is known and write an equation for it just like all the other variables, then I replace that equation by the equation \$x_{50}=1\$.

Second, this is a system of linear equations. Notice that on both sides of the equation we have \$a_{ij}b_{kij}\$ with \$\{i,j,k\}=\{s,t,r\}\$ on the left and \$\{i,j,k\}=\{r,t,s\}\$ on the right, and both are summed over \$t\$ which corresponds to \$j\$. So first I find the matrix \$C_{ik}=\sum_ja_{ij}b_{kij}\$, which requires quite some tweaking due to the messed up dimensions/slices. Now the equation reads

\$\displaystyle x_{s}\sum_{r=1}^R C_{sr} = \sum_{r=1}^R C_{rs}x_r \$

for each \$s\in1,2,\ldots,S\$. The sum of the left is just sum(C,2) in Matlab, and the sum on the right is the matrix multiplication of \$C^t\$ with the vector of unknowns \$\vec x\$. We can write the equation in matrix form as

\$\displaystyle \vec c\vec x-C\vec x=0 \$

where \$[\vec c]_s=\sum_r C_{sr}\$. Written another way, we have

\$\displaystyle \left(\mathrm{diag}(\vec c)-C\right)\vec x=0. \$

Now since we know \$x_{50}=1\$ just replace the last row of the matrix with \$[0,0,\ldots,0,1]\$. Then if \$M=\mathrm{diag}(\vec c)-C\$ our equation is \$M\vec x=\vec v\$ where \$\vec v=[0,0,\ldots,0,1]^t\$. We can solve this system directly using \.

%% Define parameters
T = 3; % Number of types (T)
S = 3; % Number of senders (S)
R = S; % Number of receivers (R), equals the number of senders
A = [.1 .1 .1   % Columns are senders (S), rows are types (T)
     .3 .3 .6
     .6 .6 .3]; 
B1 = [.9 .0 .1     
      .5 .1 .2     
      .8 .2 .1];   
B2 = [.0 .7 .2     
      .2 .9 .0     
      .1 .8 .1];   
B3 = [.1 .3 .7     
      .3 .0 .8     
      .1 .0 .8];   

C = squeeze(sum(A.*cat(3,B1,B2,B3),1)).';

M = diag(sum(C,1))-C; % matrix equations
M(end,:) = [zeros(1,S-1) 1]; % replace the last equation

v = [zeros(S-1,1);1]; % modify final RHS entry

F = M\v;

An additional benefit is that solving linear systems is much faster than solving nonlinear systems, so this code is much faster. I couldn't test with different R, S, or T so something might have to change, I'm not sure if it will just work, but I think it should.

---

With the new method for generating the equations, we can make a few more simplifications to avoid having to concatenate the Bi's explicitly, using a cell array:

T = 3; % Number of types (T)
S = 3; % Number of senders (S)
R = S; % Number of receivers (R), equals the number of senders
A = [.1 .1 .1   % Columns are senders (S), rows are types (T)
     .3 .3 .6
     .6 .6 .3]; 
B{1} = [.9 .0 .1     
      .5 .1 .2     
      .8 .2 .1];   
B{2} = [.0 .7 .2     
      .2 .9 .0     
      .1 .8 .1];   
B{3} = [.1 .3 .7     
      .3 .0 .8     
      .1 .0 .8];   

C = squeeze(sum(A.*cat(3,B{:}),1)).';

M = diag(sum(C,1))-C;
M(end,:) = [zeros(1,S-1) 1];

v = [zeros(S-1,1);1];

F = M\v;