# Raising elementwise a large symmetric, positive definite matrix to a power in Matlab

I have created a large positive definite (thus symmetric) matrix in Matlab. For sake of simplicity, let's assume that the matrix has the following form:

$K= \left[ \begin{array}{ccc} k_{11} & \cdots & k_{1n} \\ \vdots & \ddots & \vdots \\ k_{n1} & \cdots & k_{nn}\end{array} \right]$.

Now I want to raise it to a power in an elementwise way, let's say to the power of $\lambda$.

Since $K$ is known, while $\lambda$ is unknown and has to be estimated via an optimization routine, I save beforehand the matrix $K$. In the optimization routine I have to compute at each iteration the value of the matrix raised to the power of $\lambda$.

My question is the following: Is there a way to perform this operation in a fast way, by exploiting the properties of the matrix $K$?

What I have done until now:

1st approach: Compute $K^\lambda$ in a straightforward way, that is to say, using command "K.^\lambda" in matlab.

2nd approach: Use two for loops, where the first one iterates from i=1:n and the second one from j=i:n. This will create an upper diagonal matrix and due to symmetry, I can use triu command to get back the full matrix.

Nonetheless, both approaches are rather time consuming, due to the large size of the matrix. Does anyone have any better idea about how I could perform this operation?

• One thing you can do is use MATLAB's tools to time and/or profile both approaches so that you can determine empirically which approach is faster. My guess is that the K.^lambda approach will be faster, because it is vectorized. Feb 9 '15 at 23:20
• Once you have computed $K^{\lambda}$, what are you doing with it? Chances are that you can get a more efficient algorithm by restructuring the whole computation, not just this step. Mar 11 '15 at 16:02

I don't think you'll be able to do better than your second approach, but keep in mind that Matlab is column-major (like fortran):

for j=1:n
K_lam(1:j,j) = K(1:j,j).^lambda
end
K_lam = K_lam+triu(K_lam,1)'

• maybe vectorize your lower triangle beforehand then use built-in .^ can be a bit faster? Feb 9 '15 at 21:32
• +1: Matlab for loops are nowadays faster that they used to be, so this is a good approach. My own answer differs in the fact that I think that it should not be necessary to reconstruct the full matrix, but a partial storage could be sufficient for most purposes. Jun 9 '15 at 20:34

I see three possible storage strategies:

>> n = 1000;
>> l = 2.34;
>> K = randn(n,n); K = K*K';
>> Ku = triu(K);
>> p = 0; for i = 1:n, kvu(p+(1:i)) = K(1:i,i); p = p + i; end


where

• K is the full matrix
• Ku is the upper triangular part stored as a full matrix
• kvu is a vectorization of the non zero elements of Ku

and these are my timings:

>> tic, r = K.^l; toc
Elapsed time is 0.033273 seconds.
>> tic, r = Ku.^l; toc
Elapsed time is 0.029381 seconds.
>> tic, r = kvu.^l; toc
Elapsed time is 0.019505 seconds.


The Ku approach is interesting, since a lot of matrix functions for symmetric positive definite matrices access only the upper triangular part: e.g.

>> isequal(chol(K), chol(Ku))

ans =

1


Of course the ideal situation is if the optimisation algorithm could be rewritten in terms of kvu instead of K.

• Strangely enough, I answered this question thinking it was a new one, but my answer is almost 4 months late... I should double check before answering. Jun 9 '15 at 20:46
• Very thorough! Changing the storage is always fastest, but its not always worth it. +1 for quantifying the trade-off. How do the speed-ups depend on n? Jun 10 '15 at 1:10
• @MaxHutchinson: with all the data formats the problem remains $O(n^2)$, but I do expect that the "speed-up" depends also from $n$, due to a lot of factors, like memory, cache, type of results (real or complex) whose analysis is out of the scope of my answer. My intent here was only to point out that some choices of data structure (typical when linking linear algebra routines in Fortran or C) apply also for matlab scripts. Jun 11 '15 at 13:21