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

Question

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?

Answer 1

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)'

Answer 2

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.