 function Q = Rank1SeqProd(V)
 % V is an m-by-p lower triangular matrix with p<m
 % Q = V_1*V_2*...*V_p where V_k = eye(m,m) - V(:,k)*V(:,k)'
 
 [m,p] = size(V);
 
 % Method 1. Here is an explicit  follow-the-formula solution...
 Q = eye(m,m);
 for k=1:p
     v = V(:,k);
     Q = Q*(eye(m,m) - v*v');
 end
 % Flops = p*(2m^3) = 2pm^3
 Q = Q
 
 % Method 2. Here is a solution that exploits the I-plus-rank-1
 % structure of the factor matrices...
 Q = eye(m,m);
 for k=1:p
     v = V(:,k);
     z = Q*v;
     Q = Q - z*v'; 
 end
 % Flops = p*m^2*(2 + 2) = 4pm^2
 Q = Q
 
 % Method 3. Here is a solution that exploits the I-plus-rank-1
 % structure of the factor matrices and the fact that
 % V is lower triangular
 Q = eye(m,m);
 for k=1:p
     v = V(k:m,k);
     z = Q(:,k:m)*v;
     Q(:,k:m) = Q(:,k:m) - z*v'; 
 end
 % Flops = sum from k=1 to p of 4(m-k)m = 4(p*m^2 - p^2*m/2)
 Q = Q
 
 % Method 4. Here is a solution that exploits the I-plus-rank-1
 % structure of the factor matrices and the fact that
 % V is lower triangular and that it is better to
 % multiply the factor matrices in reverse order. 
 Q = eye(m,m);
 for k=p:-1:1
     % If Q [ Q11 Q12; Q21 Q22] with Q11 (k-1)-by-(k-1), then 
     % Q12 = 0, Q21 = 0 and Q11 = I. Q22 need only be involved
     % in the update...
     v = V(k:m,k);
     w = (v'*Q(k:m,k:m))';
     Q(k:m,k:m) = Q(k:m,k:m) - v*w';
 end
 % Flops = sum from k=1 to p of 4(m-k)^2 = 4(p*m^2 - p^2*m + p^3/3) 
 Q = Q
 
 % Summary if p = m...
 %    Method 1 = 2m^4
 %    Method 2 = 4m^3
 %    Method 3 = 2m^3
 %    Method 4 = (4/3)m^3
 
 % Roughly 3 points for exploiting the rank-1's and 2 points for the zeros