2
\$\begingroup\$

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}

\$\endgroup\$
  • \$\begingroup\$ 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? \$\endgroup\$ – David Jan 28 at 23:19
  • \$\begingroup\$ 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. \$\endgroup\$ – Aristide Herve Jan 29 at 7:15
4
\$\begingroup\$

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;
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.