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: 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}

• David: I've spent a bit of time on this, but your inconsistent row/column indexing is making things confusing. For A you have a_{ij} for column i row j and for B you have b_{ijk} for row i column j slice k. Is this correct, a_{11}+a_{12}+a_{13}=1? And is x_50=1?
• Aristide Herve: Yes, sum(A) is a vector of ones, as is sum(B,2). Yes, x_50 = 1 (in the small example shown in the question, x_3 = 1). Sorry for the inconsistent indexing. I added an illustration of the B array.

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;