Trefethen & Bau & MATLAB & Julia, Lecture 8: Gram-Schmidt

computing
math
teaching
Author

Toby Driscoll

Published

September 19, 2016

This lecture is about the modified Gram-Schmidt method and flop counting. The notebooks are here.

I’m lost.

Almost as an afterthought I decided to add a demonstration of the timing of Gram-Schmidt compared to the asymptotic O(n^2) flop count. Both MATLAB and Julia got very close to the trend as n got into the hundreds, using vectorized code:

n_ = collect(50:50:500);
time_ = zeros(size(n_));
for k = 1:length(n_)
    n = n_[k];
    A = rand(1200,n);
    Q = zeros(1200,n);  R = zeros(600,600); 
    
    tic();
    R[1,1] = norm(A[:,1]);
    Q[:,1] = A[:,1]/R[1,1];
    for j = 2:n
        R[1:j-1,j] = Q[:,1:j-1]'*A[:,j];
        v = A[:,j] - Q[:,1:j-1]*R[1:j-1,j];
        R[j,j] = norm(v);
        Q[:,j] = v/R[j,j];
    end
    time_[k] = toc();
end

using PyPlot
loglog(n_,time_,"-o",n_,(n_/500).^2,"--")
xlabel("n"), ylabel("elapsed time")

I noticed that while the timings were similar, Julia lagged MATLAB just a bit. I decided this would be a great chance for me to see Julia’s prowess with speedy loops firsthand.

Compare the vectorized and unvectorized Julia versions here:

Look at the last line–it’s allocating 1.4GB of memory to make the nested loop version happen! I thought perhaps I should use copy to create v in each pass, but that change didn’t help. I even tried writing my own loop for computing the dot product, to no avail.

It did help a little to replace the line in which v is updated with

v = broadcast!(-,v,Q[:,i]*R[i,j])

The bang on the name of the function makes it operate in-place, overwriting the current storage. Apparently Julia will create some syntactic sugar for this maneuver in version 0.5. Here it reduced the memory usage to 1.1 GB.

Julia’s reputation is that it’s great with loops, especially compared to MATLAB and Python. As a Julia newbie I recognize that there may still be only a small change I need to make in order to see this for myself. But I feel as though having to use that broadcast!, or even the more natural .= that may be coming, is already too much to ask. I’m frustrated, confused, and disappointed.