Set of linear equations

Question

I have \$ N^2 \$ matrices. Each one is a \$ 3 \times 3 \$ matrix. I want to concatenation them to a \$ 3N \times 3N \$ matrix. I know that the big matrix is symmetric. Evaluation of whole \$ 3 \times 3 \$ matrices are time consuming so I want to speed up my program.

Do you have any suggestion to speed it up?

clc 
load rv
load a %# Nx3 matrix [x1 y1 z1;x2 y2 z2;... ].
 %# vectors construct a sphere
L=(301:500)*1e-9; K=2*pi./L; %# 1x200 array 
%# some constants =========================================================
I=eye(3);
e0=1;
[npoints,ndims]=size(rv); 
N=npoints;
d0=(4*pi/(3*N))^(1/3)*5e-9;
V=N*d0^3; aeq=(3*V/(4*pi))^(1/3); 
E0y=ones(N,1);
E0z=E0y;
Cext=zeros(1,length(L));
Qext=zeros(1,length(L));
A=zeros(3,3,N^2);
% =================================
r=sqrt(sum(rv,2).^2); %# [Nx1] lengrh of each rv vector
for i=1:N
 for j=1:N
 dx(i,j)=rv(i,1)-rv(j,1); %# The x component of distances between each 2 point [NxN]
 dy(i,j)=rv(i,2)-rv(j,2); %# The y component of distances between each 2 point [NxN]
 dz(i,j)=rv(i,3)-rv(j,3); %# The z component of distances between each 2 point [NxN]
 end
end
dv=cat(3,dx,dy,dz); %# d is the distance between each 2 point (a 3D matrix) [NxNx3]
d=sqrt(dx.^2+dy.^2+dz.^2); %# length of each dv vector
nx=dv(:,:,1)./d; ny=dv(:,:,2)./d; nz=dv(:,:,3)./d;
n=cat(3,nx,ny,nz); %# Unit vectors for points that construct my sphere
for s=1:length(L)
 E0x=exp(1i*K(s)*rv(:,1))'; % ' #closing the quote for syntax highlighting
 % 1x200 array in direction of x(in Cartesian coordinate system)
 % Main Loop =================================================
 p=1; 
 for i=1:N 
 for j=1:N 
 if i==j 
 A(:,:,p)=a(s)*eye(3); %# 3x3 , a is a 1x200 constant array 
 p=p+1; %# p is only a counter 
 else 
 A(:,:,p)=-exp(1i*K(s)*d(i,j))/d(i,j)*(-K(s)^2*([nx(i,j);ny(i,j);nz(i,j)]... 
 *[nx(i,j) ny(i,j) nz(i,j)]-I)+(1/d(i,j)^2-1i*K(s)/d(i,j))... 
 *(3*[nx(i,j);ny(i,j);nz(i,j)]*[nx(i,j) ny(i,j) nz(i,j)]-I)); 
 p=p+1; 
 end 
 end 
 end 
%===============================================================
B = reshape(permute(reshape(A,3,3*N,[]),[2 1 3]),3*N,[]).'; %# From :gnovice
%# concatenation of N^2 3x3 matrixes into a 3Nx3N matrix
 for i=1:N
 E00(:,i)=[E0x(i) E0y(i) E0z(i)]';
 end
 b=reshape(E00,3*N,1);
 P=inv(B)*b;
 Cext(s)=(4*pi*K(s))*imag(b'*P);
 Qext(s)=Cext(s)/(pi*aeq^2);
 end
Qmax=max(Qext); Qext=Qext/Qmax;
L=L*1e9;
plot(L,Qext,'--');figure(gcf)
%# The B matrix is symmetric(B_ij=B_ji) so I can reduce the computations
%# I think I should write sth like this
% for i=1:N 
% for j=i:N 
% if i==j 
% A(:,:,p)=a(s)*eye(3); %# 3x3 , a is a 1x200 constant array 
% p=p+1; %# p is only a counter 
% else 
% A(:,:,p)=... [just same as above] 
% p=p+1; 
% end 
% end 
% end 
% But how concatenate them like befor in order?

Pls get a.mat : http://dl.dropbox.com/u/21031944/Stack/a.mat

 rv.mat: http://dl.dropbox.com/u/21031944/Stack/rv.mat

Or sth like this (for symmetry part):

for i=1:N 
 for j=1:N 
 if i==j 
 A(:,:,p)=a(s)*eye(3); %# 3x3 , a is a 1x200 constant array 
 p=p+1; %# p is only a counter 
 elseif i>j 
 A(:,:,p)=... [just same as above] 
 p=p+1; 
 else
 A(:,:,p)=A(:,:,??);
 p=p+1;
 end 
 end 
 end 

I don't know is program clear?

• Abolfazl: A bigger rv
• Abolfazl: i.sstatic.net/AJKJl.png i.sstatic.net/AJKJl.png i.sstatic.net/SOpya.png i.imgur.com/podvp.png
• Dennis Jaheruddin: Could you run the profiler to find out which part exactly requires the speedup?

Answer 1

In my experience with finite element problems, I found that creating the large matrix A easily took as long, if not longer, than calculating the solution x = A/b. Matlab does a good job with that matrix inverse.

In (old) FORTRAN creation of A is done with nested loops. In Matlab this iteration is very slow. One trick is to reverse the order of the nesting, so you iterate on small dimensions (the 3x3), and use matrix operations on the large ones (N). A further step is to do everything with matrix operations (on 4D arrays if needed), but that is harder to conceptualize.

I suspect your A matrix is sparse (lots of eye components). In that case the Matlab sparse matrix operations can be your friend. Basically you create 3 large vectors, I, J, V, where 2 identify the coordinates, and the third the (nonzero) value at each. If a specific pair of coordinates is repeated, the corresponding values are added, which in some matrix definitions is convenient. But that is more than I can explain here.

https://codereview.stackexchange.com/a/30101/27783 takes a simple heat diffusion model, using loops and indexing, and speeds it up in several steps. Simply vectorizing the inner loops helps a lot. The best speed came generating a matrix using eye (or diag), and doing fast matrix multiplication (with np.dot).

• Dennis Jaheruddin: If the goal is just solving linear equations, then indeed using ldivide or rdivide seems to be the way to go.