# matrix that forms an orthogonal basis with a given vector

A linear algebra question;

Given a k-variate normed vector u (i.e. u : ||u||_2=1) how do you construct \Gamma_u, any arbitrary k*(k-1) matrix of unit vectors such that (u,\Gamma_u) forms an orthogonal basis ?

I mean: from a computationnal stand point of view: what algorithm do you use to construct such matrices ?

The naive approach would be to apply Gram Schmidt orthogonalisation of u_0, and k-1 randomly generated vectors. If at some point the GS algorithm generates a zero vector, then you have a linear dependency in which case choose the vector randomly again.

However this method is unstable, small numerical errors in the representation of the vectors gets magnified. However there exists a stable modification of this algorithm:

Let a_1 = u, a_2,...a_k be randomly chosen vectors

```for i = 1 to k do
vi = ai
end for

for i = 1 to k do
rii = |vi|
qi = vi/rii
for j = i + 1 to k do
rij =<qi,vj>
vj =vj −rij*qi
end for
end for
```

The resulting vectors v1,...vk will be the columns of your matrix, with v1 = u. If at some point vj becomes zero choose a new vector aj and start again. Note that the probability of this happening is negligible if the vectors a2,..,ak are chosen randomly.

You can use Householder matrices to do this. See for example http://en.wikipedia.org/wiki/Householder_reflection and http://en.wikipedia.org/wiki/QR_decomposition

One can find a Householder matrix Q so that Q*u = e_1 (where e_k is the vector that's all 0s apart from a 1 in the k-th place) Then if f_k = Q*e_k, the f_k form an orthogonal basis and f_1 = u. (Since Q*Q = I, and Q is orthogonal.)

All this talk of matrices might make it seem that the routine would be expensive, but this is not so. For example this C function, given a vector of length 1 returns an array with the required basis in column order, ie the j'th component of the i'th vector is held in b[j+dim*i]

```   double*  make_basis( int dim, const double* v)
{
double* B = calloc( dim*dim, sizeof * B);
double* h = calloc( dim, sizeof *h);
double  f, s, d;
int i, j;

/* compute Householder vector and factor */
memcpy( h, v, dim*sizeof *h);
s = ( v[0] > 0.0) ? 1.0 : -1.0;
h[0] += s;
f = s/(s+v[0]);

/* compute basis */
memcpy( B, v, dim * sizeof *v); /* first one is v */
/* others by applying Householder matrix */
for( i=1; i<dim; ++i)
{   d = f*h[i];
for( j=0; j<dim; ++j)
{   B[dim*i+j] = (i==j) - d*h[j];
}
}
free( h);
return B;
}
```