title: MATLAB Optimisation of Weighted Gram-Schmidt Orthogonalisation
body: <p>I have a function in MATLAB which performs the <a href="https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process" rel="nofollow noreferrer">Gram-Schmidt Orthogonalisation</a> 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?</p>
  
  <p>I have tried converting to a MEX file but I lose parallelisation with the compiler I'm using and so it is then slower.</p>
  
  <p>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)</p>
  
  <p>Can anyone vectorise this code or make it faster? I am not sure how to do it elegantly ...</p>
  
  <p>I know the stackoverflow minds here are amazing, consider this a challenge :)</p>
  
  <h2><strong>Function</strong></h2>
  
  <pre><code>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
  </code></pre>
  
  <p>where <code>A</code> is an <code>m x n</code> matrix of complex numbers and <code>w</code> is an <code>m x 1</code> vector of real numbers.</p>
  
  <h2><strong>Bottle-neck</strong></h2>
  
  <p>This is the expression for <code>R(i,j)</code> which is the slowest part of the function (not 100% sure if the notation is correct):
  <img src="https://i.stack.imgur.com/s0tHT.jpg" alt="Weighted Inner-Product"></p>
  
  <p>where <code>w</code> is a non-negative weight function.
  The weighted inner-product is mentioned on several Wikipedia pages, <a href="http://en.wikipedia.org/wiki/Weight_function#Inner_product" rel="nofollow noreferrer">this is one on the weight function</a> and <a href="https://en.wikipedia.org/wiki/Orthogonality#Orthogonal_functions" rel="nofollow noreferrer">this is one on orthogonal functions</a>.</p>
  
  <h2><strong>Reproduction</strong></h2>
  
  <p>You can produce results using the following script:</p>
  
  <pre><code>A = complex( rand(360000,100), rand(360000,100));
  w = rand(360000, 1);
  [Q, R] = Gram_Schmidt(A, w);
  </code></pre>
  
  <p>where <code>A</code> and <code>w</code> are the inputs. </p>
  
  <h2><strong>Speed and Computation</strong></h2>
  
  <p>If you use the above script you will get profiler results synonymous to the following:
  <img src="https://i.stack.imgur.com/qVszN.png" alt="Profiler Times">
  <img src="https://i.stack.imgur.com/x4Q7l.png" alt="Profiler Code Times"></p>
  
  <h2><strong>Testing Result</strong></h2>
  
  <p>You can test the results by comparing a function with the one above using the following script:</p>
  
  <pre><code>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 );
  </code></pre>
  
  <p>where <code>Gram_Schmidt</code> is the function described earlier and <code>Gram_Schmidt2</code> is an alternative function. The results <code>zeros1</code> and <code>zeros2</code> should then be very close to zero.</p>
  
  <h2><strong>Note:</strong></h2>
  
  <p>I tried speeding up the calculation of <code>R(i,j)</code> with the following but to no avail ...</p>
  
  <pre><code>R(i,j) = ( w' * (   v    .* conj( Q(:,i) ) ) ) / ...
           ( w' * ( Q(:,i) .* conj( Q(:,i) ) ) );
  </code></pre>