# Algorithm to find an orthogonal basis (orthogonal to a given vector)

Let $K$ be a given integer, with $K$ even (and "large"). Let $\mathbf{v} \in \mathbb{R}^{K \times 1}$ be a given non-zero (column) vector. Write a (possibly efficient) algorithm to construct a matrix $\mathbf{B} \in \mathbb{R}^{K \times K-1}$ such as that: $(1)$ the column-vectors of $\mathbf{B}$ are orthogonal to each other; $(2)$ $\mathbf{B^T v} = \mathbf{0}$, where $^T$ denotes transpose and $\mathbf{0}$ a ((K-1)-dimensional) vector of all zeros.

Denoting with $\mathbf{b}_i$, $i=1,\dots,K-1$, the $i$-th column of $\mathbf{B}$, the first requirement can be rewritten as: $\mathbf{b}_i^T \mathbf{b}_j = 0$ for $i \neq j$.

Essentially, the problem asks to write an algorithm to find an orthogonal basis to the orthogonal complement of the (mono-dimensional) space spanned by $\mathbf{v}$. The algorithm can be written in pseudo-code (or in MATLAB-like or in C-like code). The algorithm takes in input the vector $\mathbf{v}$ and the integer $K$; its output is the matrix $\mathbf{B}$.

Note: the algorithm cannot do a "random trial-and-error", i.e., generate a random vector, try if it is orthogonal to $\mathbf{v}$ and to the previously found columns of $\mathbf{B}$, discard if it not, memorize it if it is. This is an explicitly forbidden brute-force approach. However, it is indeed allowed to do this as an initializing stage, i.e., for the first column of the matrix or as a "random guess" at the start, if any need to do so should arise. "Normality", i.e., finding an orthonormal basis, is not required. Input checking (e.g. checking if $K$ is an even integer) is not required.

EDIT: My thoughts and previous attempts: The problem is essentially an implementation of Gram-Schmidts orthonormalization process. However, it cannot simply be used as stated, because Gram-Schimdts' theorem assumes to start with a basis (which we do not have). What we can actually construct is a spanning set of vectors, i.e. $\mathbf{v}, e_1, \dots, e_K$, where $e_i$ denotes the $i$-th canonical base vector.

I have already implemented Gram-Schmidts, paying attention to numerical issues. The problem is a slight generalization of the process, in which you do not start with a basis, but with a spanning set, find a suitable basis in the set (I don't know how to do it), which must contain $\mathbf{v}$ as its first element, and then apply Gram-Schimdts process.

P.S. Any help is appreciated, I am no master of linear algebra, especially numerical implementations. Thanks.

• The onus is on you to test, you get what you pay for. Mar 12, 2014 at 23:06
• For orthogonalizing $v,e_1,\ldots,e_K$, you would need to implement some sort of rank revealing Gram-Schmidt (see en.wikipedia.org/wiki/RRQR_factorization); SVD can be useful or simply Gram-Schmidt with pivoting could work as well. I just wonder how that "random trial-and-error" approach would work when one normally does not have the eternity available for calculations. Mar 13, 2014 at 12:00
• @AlgebraicPavel: I used to have a 'fear' of random algorithms (I'm thinking continuous parameter optimization here), but having seen the amount of randomness added by spurious programming/design (and things still work), I have moved beyond my fear :-). Mar 13, 2014 at 16:07
• @AlgebraicPavel: It is exactly the reason why it was explicitly forbidden (although it was obvious). Note, however, that for extremely low $K$, it does work and that, for very big $K$, a "rationalized" random trial-and-error (with optimization) could actually work, since, for $K$ big, vectors tend to naturally be "numerically orthogonal" to each other. With optimization, I underline. "And things still work", strangely enough. :-) Mar 13, 2014 at 20:45

By the way, most likely a much simpler and efficient way to compute your matrix $B$ here is just using a single Householder reflection. Consider $w=v\pm\|v\|_2e_1$ (the sign in "$\pm$" is usually chosen to be same as in the first component in $v$ from numerical reasons). Then set $Q=I-2ww^T/w^Tw$. Since $Qv=\mp\|v\|_2e_1$ (depending on the choice of the sign in $w$), we have $v=\mp\|v\|_2Qe_1$ (multiply by $Q^T=Q$). It means that $v$ is a scalar multiple of the first column (row) of the orthogonal matrix $Q$. Let $Q=[q_1,\ldots,q_K]$. It follows that $q_i^Tv=0$ for $i=2,\ldots,K$. Hence $B=[q_2,\ldots,q_K]$ satisfies $B^Tv=0$ and $B^TB=I$.

In contrast to the $O(K^3)$ complexity of the Gram-Schmidt approach (which is of course generally applicable to compute orthonormal basis for any given set of vectors), this way you get your $B$ using $O(K^2)$ operations if you need the full matrix explicitly computed or $O(K)$ if it's enough to have it implicitly represented by $w$.

• A lot of time has passed, but your answer deserves compliments. It is, in fact, an efficient algorithm and the one I'm actually using is (a variant of) this one. Initially, I was unaware of "Householder magic" (the unofficial name of Householder transformations), thus I went with QR/MGS and was satisfied with that. Thank you for your contribution. Mar 27, 2014 at 22:06
• This is a beautiful and elegant solution. Did you get this from some book? Or you thought of this yourself? Dec 11, 2020 at 14:40
• There's a question on meta asking to clarify this answer : math.meta.stackexchange.com/questions/32971/… Jan 3, 2021 at 12:41
• @TeresaLisbon: Interesting. I’d never imagine that :-) I don’t read meta very often. Jan 3, 2021 at 13:24
• @NadaBrothers: Thanks, credits go to Householder :-) Jan 3, 2021 at 13:26

Run Gram Schmidt orthogonalization on the vectors $v,e_1,...,e_K$ (ignore one of the reductions as it will be zero or near zero if numerics are an issue), to get $v_1,...,v_K$. Then let $B= \begin{bmatrix} v_2 \cdots v_K \end{bmatrix}$.

Details:

By $e_k$ I mean the vector of zeroes with a one in the $k$th position. Note that $e_1,...,e_K$ forms a basis for $\mathbb{R}^K$.

Let $u_1 = v, u_k = e_{k+1}$ for $k=1,...,K$. Note that $u_1,...,u_{K+1}$ is a spanning set but is not a basis.

Then the algorithm (this is theoretical, in a practical situation the comparison in Step 3 would be replaced by something like $\|w_i\| < \epsilon$) looks like:

1. $v_1 = {u_1 \over \|u_1\|}$, $i=2$
2. $w_i = u_i - \sum_{j=1}^{i-1} v_j^T u_i$
3. if $w_i \neq 0$ then $v_i = { w_i \over \|w_i\|}$ else $v_i = 0$
4. $i \leftarrow i+1$, if $i \le K+1$ goto Step 2

This is basically the Gram Schmidt algorithm, except one of the $v_i$s will be zero.

By construction, $\operatorname{sp} \{ u_1,...,u_i \} = \operatorname{sp} \{ v_1,...,v_i \}$, hence the final set will span $\mathbb{R}^K$ (since $u_2,...,u_{K+1}$ is a basis). Consequently exactly one of the $v_i$s will be zero.

Also by construction, $v_j^T v_i = 0$ for all $j<i$, and $\|v_i\| = 1$ for all $i$ except one. Let $v'_i, i=1,...,K$ be the non-zero $v_i$s.

Now let $B = \begin{bmatrix} v'_2 \cdots v'_K \end{bmatrix}$ and note that $B^T v'_1 = 0$, and since $v = \|v\| v'_1$, we are finished.

• Thanks! It (partially) answers my question: I see a problem in that. Gram-Schmidts orthogonalization assumes you already have a basis to start with (I recently checked it: the theorem starts with a basis, not a spanning set). If, for example, $\mathbf{v}$ = $\mathbf{e}_1$, then what is the "discard criterion" when you do the scalar products $\mathbf{v}^T \mathbf{e}_i$? P.S. When I said "test" I didn't mean numerical tests (of which the onus is on me), I apologize for the poor wording. What I actually meant was simply theoretically correct algorithms. Mar 13, 2014 at 8:39
• @PseudoRandom: I have added details above. There is no problem, the algorithm basically discards vectors that are in the span of the previously computed elements. This is a theoretical algorithm, it would need to be adapted to handle numerics (or better still, just use a QR decomposition). Actually, I am always a little surprised that the Gram Schmidt algorithm and minor extensions as above are not emphasised more in classes. It is a very useful tool numerical or not. Mar 13, 2014 at 16:01
• In other words, you just proved that, in the Gram-Schimdt theorem, the hypothesis of starting with a basis is uncalled for: a spanning set is more than enough. This was unknown to me. Thank you! P.S. I personally don't like QR decomposition. Numerically speaking, it is not as good as it seems (I am aware of MGS, Modified Gram Schmidt). If you only have to find an orthogonal basis, the normalizing step is not a great idea. It often occured to me that orthogonality (even approximate orthogonality) is lost with normalization, but may be keeped if you do not normalize. Anyway, thanks! Mar 13, 2014 at 20:57
• Well, the usual Gram Schmidt presentation doesn't account for the case when the initial collection is linearly dependent, but (numerics aside) this is easy to deal with. The modified method (MGS, that is) is generally taken to have better properties, the above can be easily modified to avoid the premature normalizing. Mar 13, 2014 at 21:07
• Question: can the normalization step simply be skipped altogether, since the goal is an orthogonal basis (and not an orthonormal basis)? If (as I think), yes, it could help orthogonality, numerically speaking. P.S. We never considered any variation (no matter how minor) of the Gram-Schimdt process. Mar 13, 2014 at 21:14