Posted: April 18, 2019
Instead of re-inventing the wheel, I googled the code for the Gram-Schmidt process. The first result is a note from MIT, on how to code the G-S process in MATLAB.
The code is pretty straight forward.
n column vectors.
Q stores the normalized orthogonal vectors, and
R is the $n\times n$ matrix such that $A=Q\cdot R$.
for j=1:n v=A(:,j); for i=1:j-1 R(i,j)=Q(:,i)'*A(:,j); % change A(:,j) to v for more accuracy v=v-R(i,j)*Q(:,i); end R(j,j)=norm(v); Q(:,j)=v/R(j,j); end
But when I coded it in Fortran, weird things start to happen.
do j=1,7 vec=Amat(:,j) do i=1,j-1 Rmat(i,j) = dot_product(Qmat(:,i),vec(:)) vec(:)=vec(:)-Rmat(i,j)*Qmat(:,i) ! subtract the projection end do ! now vec the perpendicular to all Q(:,i) Rmat(j,j) = cdsqrt(dot_product(vec(:),vec(:))) Qmat(:,j)= vec(:)/Rmat(j,j) ! normalize end do
The above code does not give the expected result.
Moreover, when I check the rank of the
Amat, it gives 3. But after orthogonalization, the rank of
Qmat is 7!
call zgesvd("N", "N", dm, 7, Amat, dm, sing, Umat, dm, vt, 7, Awork, 25*dm, Arwork, inf) write(*,*) "************************Rank****************************" rank = 0 do i=1,7 if (sing(i).ge. 0.05d0) rank = rank +1 end do write(*,*) inf
This is when I started to question everything. Is it that the
zgesvd not for calculating the rank? Or is it that the process is not as such? I checked and checked everything again and again, making unit tests, giving it dummy input and everything seems fine. But the inconsistency keeps arising.
I finally realized that the problem arises when the input are “realistic”, i.e., with limited precisions. When the inputs are “nice” such as
the results are correct anyway. But once you apply the above code to vectors such as
even if the two vectors are practically linearly dependent, the error in calculation will result in a tiny non-zero vector, which will later be normalized regardless of it’s magnitude. To avoid this, check
abs(dot_product(vec(:),vec(:)) to remove such problem.
do j=1,7 vec(:)=Amat(:,j) if (j.gt.1) then do i=1,j-1 vec(:)=vec(:)-dot_product(Qmat(:,i),vec(:))*Qmat(:,i) ! subtract the projection end do ! now vec the perpendicular to all Q(:,i) end if if (abs(dot_product(vec(:),vec(:))).le.0.01) then Qmat(:,j)=vec(:) else Qmat(:,j)= vec(:)/sqrt(dot_product(vec(:),vec(:))) ! normalize end if end do
Now the code gives the correct result.
See this thread.
I used to read the different coding style between industry-grade code and academic-project code. I have always taken pride in my code and considered the industrial way of coding ugly and ungraceful. I actually considered myself quite aware of errors by machine precisions, and have applied quite many “academic” ways to cope with errors anyway. Besides, I have never really met any problem that specifically requires me to address such error.
I fell into the pitfall mainly because
In hindsight, I should have realized the code from note from MIT is purely academic-project grade code. You can see the example in that code. The example actually applies the code to
I wouldn’t say that my version of G-S is industrial grade, but this is surely a lesson: actual programs’ behaviors are much more complicated and a pure “theoretical” way of doing it might not work. Do not ever forget the error even if someone might have done the dirty work for you.