Given a sequence of vectors $(a_1, a_2, ...)$ with $(\forall i \in \{ 1, 2, ... \} ) a_i \in \mathbb{R}^n$, we define the following sequence of matrices $(M_1, M_2, ...)$ with $(\forall i \in \{ 1, 2, ... \} ) M_i \in \mathbb{R}^{n \times n}$:
$$ \begin{aligned} M_1 &= I \\ M_{i+1} &= M_i - \frac{M_i a_i a_i^T M_i}{a_i^T M_i a_i} \\ \end{aligned} $$
Looking at the first few elements, we have:
$$ \begin{aligned} M_1 &= I \\ M_{2} &= I - \frac{a_1 a_1^T}{a_1^T a_1} \\ M_{3} &= I - \frac{(a_1^T a_1) a_2 a_2^T - (a_1^T a_2) a_2 a_1^T - (a_1^T a_2) a_1 a_2^T + (a_2^T a_2) a_1 a_1^T}{(a_1^T a_1)(a_2^T a_2) - (a_1^T a_2)(a_2^T a_1)} \\ \end{aligned} $$
The problem is to find a closed form solution for $M_i$ in terms of $(a_1, a_2, ..., a_i)$.
Note that it easy to prove by induction that $(\forall i \in \{ 1, 2, ... \} ) M_i^T=M_i$ and $M_i M_i = M_i$.
I've made the following conjecture about the solution to the problem:
$$ M_{i} = I - \frac{ \sum_{ \sigma \in S_{i-1} }{ \left[ \operatorname{sgn} (\sigma) \prod_{l=1}^{i-1}{ \left[ a_l^T a_{\sigma (l)} \right ] } \sum_{k=1}^{i-1}{ \left[ \frac{a_k a_{\sigma (k)}^T}{a_k^T a_{\sigma (k)}} \right] } \right] } }{ \sum_{ \sigma \in S_{i-1} }{ \left[ \operatorname{sgn} (\sigma) \prod_{l=1}^{i-1}{ \left[ a_l^T a_{\sigma (l)} \right ] } \right] } } $$
If my conjecture is true then it should be possible to prove it by induction, but I'm having a hard time with the algebra.
(EDIT: I'm taking a break from this problem for now, but as a note-to-self for if/when I return to it, a useful approach might be to try to isolate the $S_{i-1}$ terms from $M_{i+1}$, use this to find a simplified expression for $M_i - M_{i+1}$, and compare this expression to $\frac{M_i a_i a_i^T M_i}{a_i^T M_i a_i}$ to see if they're provably equivalent).
For context, this is part of a broader problem which is to find a closed form solution for the volume of a parallelotope defined by the vectors $(a_1, a_2, ..., a_m)$ using Gram-Schmidt orthogonalisation, which in turn is to be used to derive the Leibniz formula for the determinant of a matrix. Hence, I am looking for solutions to this problem which avoid relying on definitions or properties of determinants, as this would lead to circular reasoning in the broader problem.
In your posting, note that $M_2$ is orthogonal projection onto $(a_1/|a_1|)^\perp$
(A) We want to show that $M_{i+1}$ is orthogonal projection onto orthogonal compliment of subspace generated by $a_k, \ 1\leq k\leq i$ by induction.
So we assume that $M_i$ satisfies the assumption.
(1) Hence easily we can showw that $M_{i+1}a_k=0$ for $1\leq k\leq i$.
(2) If $v\cdot a_k=0$ for $1\leq k\leq i$, then $M_{i+1}v=v$ :
First note that $ M_iv=v$. And since $M_ia_i$ is a linear combination of $a_k,\ 1\leq k\leq i$, then $$(M_ia_i)\cdot v=0,$$ which implies $M_{i+1}v=v$. Hence we completes the proof.
(B) Now we will construct $M_{k+1}$ : Define $A = [a_1\cdots a_k]$. By construction of all $M_i$ we note that $A$ has rank $k$.
If $\pi$ is orthogonal projection onto column space of $A$, then for any $v\in \mathbb{R}^n$, there is unique $x$ s.t. $\pi v=Ax$.
Hence $$A^TA x = A^T \pi v = A^Tv $$ so that $$ x = (A^TA)^{-1} A^T v $$ which implies that $$ \pi v =Ax = A(A^TA)^{-1} A^T v $$ Hence $$ M_{k+1} = I- \pi = I- A(A^TA)^{-1} A^T $$