Create a unitary matrix out of a column vector I need an algorithmic way, to create a unitary matrix out of a first, given, column vector. i.e. if the column vector was $(0 1 1 1)^T$, I'd need an algorithm to complete it to a matrix like this.
$$
\frac{1}{\sqrt{3}}  \quad
\begin{pmatrix}
0 & 1 & 1 & 1 \\
1 & 0 & 1 & -1 \\
1 & -1 & 0 & 1 \\
1 & 1 & -1 & 0
\end{pmatrix}$$
The content of the other column vectors is not relevant, as long as it makes the matrix unitary. The only partial answer I found was in this previous post, but it only covers vectors containing only $1$.
Does anyone know a smart way to accomplish this?

EDIT
Based on the answer by @ben-grossmann, I've written this Python function to create the matrix, if anyone else needs it.
def create_unitary(v):
    dim = v.size
    # Return identity if v is a multiple of e1
    if v[0][0] and not np.any(v[0][1:]):
        return np.identity(dim)
    e1 = np.zeros(dim)
    e1[0] = 1
    w = v/np.linalg.norm(v) - e1
    return np.identity(dim) - 2*((np.dot(w.T, w))/(np.dot(w, w.T)))

 A: One solution that I find elegant is to use the Householder matrix with the correct first column. One way to build this matrix is as follows.
Let $v$ denote the desired first column. Let $e_1$ denote the column-vector $e_1 = (1\ \ 0\ \ \cdots \ \ 0)^T$. First of all, if $v$ is a multiple of $e_1$, then it suffices to either use $I$ or $-I$, where $I$ denotes the identity matrix. If not, then let $w$ denote the vector $w = v/\|v\| - e_1$. The Householder matrix corresponding to the reflection across the hyperplane orthogonal to $w$ is given by
$$
H = I - 2 \frac{ww^T}{w^Tw}.
$$
Applying the process to $v = (0,1,1,1)^T$ leads to the matrix
$$
H = \left[\begin{matrix}
0 & \frac 13 \sqrt{3} & \frac 13 \sqrt{3} & \frac 13 \sqrt{3}\\
\frac 13 \sqrt{3} & \frac 23 & -\frac 13 & -\frac 13\\
\frac 13 \sqrt{3} & -\frac 13 & \frac 23 & -\frac 13\\
\frac 13\sqrt{3} & -\frac 13 & -\frac 13 & \frac 23\end{matrix}\right].
$$

Another approach that bears mentioning: we could do the following using the null_space method.
import numpy as np
import scipy.linalg as la

def create_unitary(v):
    return np.hstack((v,la.null_space(v.T)))

A: One way to do this is to start by extending the given column to any invertible square matrix $A$ (for example by adding standard basis vectors). You can then apply the Gram–Schmidt process to the columns of $A$ to obtain an orthogonal (unitary) matrix $U$.
A: In addition to the vector of interest $(v)$, choose another arbitrary vector $(a)$. Use this vector pair to create an skew-hermitian matrix $(S)$, then use a Cayley Transform to generate the desired unitary matrix $(U)$
$$\eqalign{
S &= av^\dagger - va^\dagger \\
U &= (I+S)^{-1}(I-S)\\
}$$
The Householder transform suggested in the other answers is also a good idea.
