I am trying to vectorize a loop in which each column vector of a 2D matrix (n-by-n) is found by multiplying each single element in a diagonal matrix with a column vector in another n-by-n 2D matrix. It seems like it would be simple to vectorize, but I must be missing something. I tried converting the diagonal matrix to a column vector first, but it still needs to do the operation one element at a time.

Thank you in advance.

%% example matrix
B = [ 3     2      .9    2
      2      4       1    2
      3      4      -1    0
     .5     .5      .1    1];

%% find eigenvectors and eigenvalues
 [ve, va] = eig(B,'nobalance');

%% get size of B ... 4

%% time constants
    t = 2000e-6;

%% pre-allocate memory for Matrix Mn (Mnew) and Mp (Mprime)

%% Original code, runs correctly, not vectorized - steps through each column
for m = 1:q,
    Mn(:,m) = exp( va(m,m) * t0) * ve(:,m);             
    Mp(:,m) = exp( va(m,m) * t ) * ve(:,m);            

% quick explanation ^ for each column of Mn and Mp, a single value
% exp(va(i,i)*t) is multiplied by a column of ve(:,i)
% vectorizable?
% Mn(column) = A * ve(column), but A is an individual element in a different vector already

%% New code... vectorized?
% random starting idea:
va1 = va * [1; 1; 1; 1];  % changes diagonal matrix into column vector of values, BUT still has to calculate one at a time...

First you want the eigenvalues as a vector, i.e.

va = diag(va);

or just call eig with an output option:

[ve, va] = eig(B,'nobalance','vector');

Then to multiply the columns you can use either diagonal matrix multiplication

Mp = ve*diag(exp(va*t));

or bsxfun:

Mp = bsxfun(@times,ve,exp(va*t)');
  • $\begingroup$ Is A*diag(d) really optimized to take $O(n^2)$ operations only? I have never seen it mentioned in the documentation. $\endgroup$ – Federico Poloni Nov 11 '15 at 12:24
  • $\begingroup$ The diagonal matrix should be sparse really. For a one-off, it is probably better to do the bsxfun approach, to avoid the overhead of sparse(diag(d)). $\endgroup$ – GeoMatt22 Nov 11 '15 at 14:50

The results of the following codes are equivalent to the results of for loop

Mp=(exp( diag(va) * t ))' .* ve;
Mn=(exp( diag(va) * t0 ))' .* ve;
  • $\begingroup$ I think I must be missing something, as the operation .* does not allow for multiplying the transposed vector (exp( diag(va) * t0 )) by the 4x4 matrix ve. "Error using .* Matrix dimensions must agree." $\endgroup$ – quantumquack Nov 9 '15 at 20:46
  • $\begingroup$ Ohh sorry, I got used to Octave broadcasting. I don't have access to matlab. The above codes work properly in Octave. For matlab as @user18173 mentioned you can use bsxfun function. $\endgroup$ – Ömer Nov 10 '15 at 16:01

Your Answer

