(If you don't like the generality below, feel free to take $\mathbb{K} = \mathbb{R}$ and $b$ positive definite.)
Let $V$ be a finite-dimensional vector space over a field $\mathbb{K}$ and let $b = \langle \cdot, \cdot \rangle$ be a symmetric bilinear form on $V$; denote $q$ the associated quadratic form. Assume that $b$ is definite, i.e. $q$ has trivial isotropic cone, equivalently $b$ is nondegenerate in restriction to any subspace $W \subseteq V$. This implies that $V = W \oplus W^\perp$, in particular the orthogonal projection onto any subspace is well-defined.
There are two famous algorithms to find an orthogonal basis, starting with any basis:
- The Lagrange method (also called Gauss reduction in French, strangely): In coordinates, write $q$ as a sum of squares of linear forms by eliminating the mixed terms.
- The Gram-Schmidt process: put $e_1 = v_1$, and compute $e_{k+1}$ as the orthogonal projection of $v_{k+1}$ on the subspace spanned by $(e_1, \dots, e_k)$. Note: I do not scale the vectors to make them have unit norm.
I claim that when both methods are applicable ($b$ is definite), they are dual to each other, in the sense that Lagrange yields a change of coordinates and Gram-Schmidt a change of basis that correspond to each other. I think this is easy to see on a simple example.
I am surprised to not find this connection spelled out anywhere, especially in a book/lecture notes that covers both. But maybe I just didn't search well enough. Have you heard this before?