In Digital Signal Processing (DSP) one should often work or encounter the notion of vector and function spaces and finally orthonormal basis. The latter is one of those small building blocks which makes so many devices working and kicking.

The goal of this article is to publish MATLAB source code of the Gram-Schmidt algorithm. Given a set of linearly independent vectors, the algorithm returns a set of orthonormal basis. It is said that this basis spans the vector space. In other words, any vector from this vector space can be uniquely expressed using those basis vectors.

This algorithm enforces the so-called orthogonality principle. This principle is one of corner stones of the modern DSP. It’s role can be hardly overestimated. For example, without it FFT wouldn’t work! And without FFT… we would be somewhere in the ages before 19th century…

The shared algorithm is the classic one. A word of caution here – the classic algorithm is numerically unstable. One should implement the modified version, which requires only a small modification.

The source code can be downloaded [HERE]. I tried to comment it heavily so that it is maximally clear what is happening at each step.

In applications where the input matrix is large and requires numeric stability, and speed, one should instead use the optimized QR decomposition routines [wiki]. The QR decomposition returns the Q matrix (exactly a matrix of orthonormal basis as for the Gram-Schmid algorithm) + an upper triangular matrix with dot products (basis i with input vector j)…

Let me share a small gem. Suppose you are give 4 out of 5 orthonormal vectors that span 5 dimensional vector space. You would like to find the 5th missing orthonormal vector. With this routine you can find it easily the following way:

- Initialize the input matrix with zeros
- Put in the first 4 columns the orthonormal vectors you have
- In the last column put just random values
- Run the routine

In the output you should get a matrix with the first 4 columns identical to those 4 given orthonormal vectors, while the last 5th column will be the missing orthonormal vector – the one you were looking for. Actually you can have several missing orthonormal vectors and the algorithm will still work.