# Vectorizing equations for fsolve

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

• 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? Jan 28, 2020 at 23:19
• 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. Jan 29, 2020 at 7:15

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;