MATLAB Optimisation of Weighted Gram-Schmidt Orthogonalisation

Question

I have a function in MATLAB which performs the Gram-Schmidt Orthogonalisation with a very important weighting applied to the inner-products (I don't think MATLAB's built in function supports this). This function works well as far as I can tell, however, it is too slow on large matrices. What would be the best way to improve this?

I have tried converting to a MEX file but I lose parallelisation with the compiler I'm using and so it is then slower.

I was thinking of running it on a GPU as the element-wise multiplications are highly parallelised. (But I'd prefer the implementation to be easily portable)

Can anyone vectorise this code or make it faster? I am not sure how to do it elegantly ...

I know the stackoverflow minds here are amazing, consider this a challenge :)

Function

function [Q, R] = Gram_Schmidt(A, w)
    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));
    v = zeros(n, 1);

    for j = 1:n
        v = A(:,j);
        for i = 1:j-1
            R(i,j) = sum(   v    .* conj( Q(:,i) ) .* w ) / ...
                     sum( Q(:,i) .* conj( Q(:,i) ) .* w );
            v = v - R(i,j) * Q(:,i);
        end
        R(j,j) = norm(v);
        Q(:,j) = v / R(j,j);
    end
end

where A is an m x n matrix of complex numbers and w is an m x 1 vector of real numbers.

Bottle-neck

This is the expression for R(i,j) which is the slowest part of the function (not 100% sure if the notation is correct):

where w is a non-negative weight function. The weighted inner-product is mentioned on several Wikipedia pages, this is one on the weight function and this is one on orthogonal functions.

Reproduction

You can produce results using the following script:

A = complex( rand(360000,100), rand(360000,100));
w = rand(360000, 1);
[Q, R] = Gram_Schmidt(A, w);

where A and w are the inputs.

Speed and Computation

If you use the above script you will get profiler results synonymous to the following:

Testing Result

You can test the results by comparing a function with the one above using the following script:

A = complex( rand( 100, 10), rand( 100, 10));
w = rand( 100, 1);
[Q , R ] = Gram_Schmidt( A, w);
[Q2, R2] = Gram_Schmidt2( A, w);
zeros1 = norm( Q - Q2 );
zeros2 = norm( R - R2 );

where Gram_Schmidt is the function described earlier and Gram_Schmidt2 is an alternative function. The results zeros1 and zeros2 should then be very close to zero.

Note:

I tried speeding up the calculation of R(i,j) with the following but to no avail ...

R(i,j) = ( w' * (   v    .* conj( Q(:,i) ) ) ) / ...
         ( w' * ( Q(:,i) .* conj( Q(:,i) ) ) );