how to replace the following code, to get the matrix

Hi, there! I going to structure the following matrix, which has some regulations.please see the picture below: if n=4;

if n=5;

My current code is below, which has many 'for'loops.

if true
  % code
    function P1=transition_probabilityP1(n)
    A1={zeros(n^2,n^2)};
    for i=1:n-1
        for j=i+1
        A1{i,j}=diag(0.407*ones(1,n));
    end;
  end;

     for i=2:n
         for j=i-1
         A1{i,j}=diag(0.298*ones(1,n-1),-1);
         end
     end   

      for i=1:n
        for j=i
        e=diag(0.295*ones(1,n-1),1);
          if i==1  
           f=diag(0.298*ones(1,n));
           A1{1,1}=e+f;
           A1{1,1}(n,n)=0.593;
          elseif i==n
            g=diag(0.407*ones(1,n));
            A1{n,n}=e+g;
            A1{n,n}(1,1)=0.705;
            A1{n,n}(n,n)=0.702;
         else 
            A1{i,j}=e;
            A1{i,j}(1,1)=0.298;
            A1{i,j}(n,n)=0.295;
         end;
      end
    end
       ind=cellfun('isempty',A1);
       [r,c]=find(ind==1);
       for i=1:length(r)
          A1{r(i),c(i)}=zeros(n,n);
       end   
       P1=cell2mat(A1);
 end

Is there anyone who know how to make the same matrix with less loops, less memory and compute faster. Generally, I want compute the matrix P with dimension larger than 3,200,000 * 3,200,000 ! Thank you！

[C.J. Harris]

A 3,200,000 x 3,200,000 matrix equals 1.024 * 10^13 elements. Assuming you want to use doubles, that's 8 bytes per element, i.e. 8.192 * 10^13 bytes total, or 81.92 TB. That's not going to work on any sort of machine you can get your hands on.

Think smaller.

[Ralf]

Since the Matrix contains many zeros, consider using the sparse Matrix concepts.

See help sparse.

[Jason]

Thank you, Ralf ! I have tried to use sparse function to build the matrix. The code is bellow :

    if true
    % code
    tic;
    u1=3.742; u2=2.707;lambda=2.741;
    beta=u1+u2+lambda;
    a=u1/beta;
    b=u2/beta;
    c=lambda/beta;
    n=4;
    A1=cell((n+1),(n+1));
    for i=1:n;
    A1(i,i+1)={a*speye(n+1)};   % why cannot use vectorize?
    end; 

    for i=2:n+1;
    A1(i,i-1)={spdiags(c*ones(n,1),-1,n+1,n+1)};
    end;

    for i=1:n+1;
    e=spdiags(repmat(b*ones(n,1),2,1),1,n+1,n+1);   
    if i==1
       A1(i,i)={e+c*speye(n+1)};
       A1{i,i}(n+1,n+1)=b+c;
    elseif i==n+1
        A1(i,i)={e+a*speye(n+1)};
        A1{i,i}(1,1)=a+c;
        A1{i,i}(n+1,n+1)=b+c;
    else
        A1(i,i)={e};
        A1{i,i}(1,1)=c;
        A1{i,i}(n+1,n+1)=b;
      end;
    end;

    ind=cellfun('isempty',A1);
    [r,c]=find(ind==1);
    for i=1:length(r)
       A1(r(i),c(i))={sparse(n+1,n+1)};
    end;

    clear u1 u2 lambda a b c n e ind r c i beta
    P1=cell2mat(A1);
    toc
  end

• 1. The weird thing is this code seems use more memory than changes the matrix directly to sparse matrix. (1900bytes vs. 1300bytes).
• 2. do you know how to reduce the ''for'' loops in the current code?

Answer 1

[Brendan Hamm]

You really just need to make a sparse diagonal matrix. This can very easily be done with spdiags. Use the last syntax: A = spdiags(B,d,m,n), which you would be calling as: A = spdiags(B,d,n^2,n^2);. The issue now is only in setting up B and d which really just requires you to replace a few values in each of the columns of B, but this can be done by finding the appropriate indices to replace (note that not all values in B will be used, so start your counting from the bottom as this function will omit unnecessary values at the top of the matrix B). Put your head to the grindstone and I'm sure you can find the patterns necessary to do this for arbitrary values of n, without loops or cell arrays. I will do it for just for the lowest diagonal:

n = 5; % Change it to 4 and verify the code works in both cases.
B = 0.298*ones(n^2,1); % Use the calculation of wherever 0.298 came from if you wish
d = -(n+1); % The diagonal where the 1st vector of B should appear
B(end-n:-n:1) = 0;
A = spdiags(B,d,n^2,n^2); % This has the lower diagonal
Q = full(A); % For visualizing the correct results.
open Q

Now, just add the columns for the other diagonals to B and add the diagonal identifier to d. Also, never use the line of code:

if true
   %code
end

true is always true so the code always runs.

[Jason]

Thank you very much, Hamm! your answer help me a lot!

[Jason]

Hi, Hamm! I have used the sparse indexing expression to build a sparse matrix, but the matlab remind me this expression is likely to be slow. Do you have any other methods to solve this problem? Thank you in advance!

if true
  % code
        Aeq1=spalloc((n+1)^2+1,(n+1)^2,(n+1)^4);
        f1=spalloc((n+1)^2,1,(n+1)^2);
        for i=1:S;          
            Aeq1(:,i)=C{i}(:,1);
            f1(i)=f{i}(1);
        end;    
end

[Brendan Hamm]

This does not look like what I suggest above. If you are creating storage for a sparse matrix, the third input is the maximum number of non-zero elements you expect and therefore determines how much memory is allocated. What you are doing is allocating nearly all of the elements of the matrix to be non-zero so takes up a lot of space. Furthermore it is going to be slower than using the spdiags approach I mention above.

[Jason]

Actually, this is another part of the code, the first part I have already done using spdiags function. In this part, I want to use sparse indexing to creat another matrix,in which C and f are cell matrix and initially I want Aeq1 and f1 choose the coulmn one of each sub-cell matrix. And I know this sparse indexing method could be slow and I am thinking if there is any other method to replace it?

if true
  % code
     C=mat2cell(Aeq,size(Aeq,1),A*ones(1,S));
     f=mat2cell(f,A*ones(S,1),1);   
     Aeq1=spalloc((n+1)^2+1,(n+1)^2,(n+1)^4);
     f1=spalloc((n+1)^2,1,(n+1)^2);
       for i=1:S;          
           Aeq1(:,i)=C{i}(:,1);
           f1(i)=f{i}(1);
       end;       
end

[Jason]

And I think the preallocate is better than do nothing.

[Stephen23]

@Jason: when you click the {} Code button without selecting some text first then it creates this placeholder text:

if true
  % code
end

However this is simply placeholder text: you do not need to keep it, it is just supposed to show where the code text should go. It is not markup! The best way to format code on this forum is simply to (in this order):

1. select the code text
2. click the {} Code button

Alternatively simply put two space characters at the beginning of every code line.

http://www.mathworks.com/matlabcentral/answers/help/markup

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup