Compute the normal vector of a hyperplane in $R^n$

Well, basically any vector $N$ that holds $N \dot (p_i - p_j) = 0$ for for any vector $p_i, p_j \in H$ $H$ is the hyperplane.

Only $n$ different points are needed, call $P$ the set of $n$ points in $H$, and define $P^* = \{ p^*_i = p_i - p_{i+1} : p_i,p_{i+1} \in H , i \in [1,n-1]\}$.

$P^*$ is the set of vectors on $H$.

Any normal vector $N$ that holds that $N \dot p^*_i , i \in [1,n-1]$ is a normal vector.

$N$ has $n$ components, so we got here $n$ unknowns. And with the $n_1$ $p^*_i$ vectors we got $n-1$ equations.

An extra equation could be defined if we want a normalized vector, but this add a quadratic term in the systems, as the vector with all zero components can't be a normal vector of any hyperplane at least one of its components has a non-zero value, suppose that the $j$ component of $N$ has value $1$. and we have then a Ax = B, with n-1 equations and unknowns. we just have to try which j makes the systm non singular and solve it.


The system of equations:
$\begin{pmatrix}
p^*_{1,1} & p^*_{1,2} & \cdots & p^*_{1,j-1} & p^*_{1,j+1} & \cdots & p^*_{1,n} \\ p^*_{2,1} & p^*_{2,2} & \cdots & p^*_{2,j-1} & p^*_{2,j+1} & \cdots & p^*_{2,n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
p^*_{n,1} & p^*_{n,2} & \cdots & p^*_{n,j-1} & p^*_{n,j+1} & \cdots & p^*_{n,n} \\
\end{pmatrix} \begin{pmatrix} N_1 \\ N_2 \\ \vdots \\ N_{j-1} \\ N_{j+1} \\ \vdots \\ N_n \end{pmatrix} = - \begin{pmatrix} p^*_{1,j} \\ p^*_{2,j} \\ \vdots \\ p^*_{(j-1),j} \\ p^*_{(j-1),j} \\ \vdots \\ p^*_{n,j} \end{pmatrix}
$

Start checking for $j \in [1,n]$ and you will get at much $n$ and at least 1 normal vector.

Python:

from numpy import array, dot
from numpy.linalg import solve, norm


def normal(P):
    n = len(P)

    if n < 1:
        raise Exception, "O-size list."

    m = len(P[0])

    if array([len(p) != m for p in P]).any():
        raise Exception, "Inconsistent set of points."

    if n != m:
        raise Exception

    P_ = [p-q for p,q in zip(P[:-1],P[1:])]

    N = []

    for i in range(n):
        try:
            x = solve(
                array([list(p)[:i] + list(p)[i+1:] for p in P_]),
                array([-p[i] for p in P_])
                )

            x = array(list(x)[:i] + [1] + list(x)[i:])
            N.append(x)

        except:
            pass

    n = N[array([(dot(n,n)-1)**2 for n in N]).argmin()]

    return n / norm(n)



if __name__ == '__main__':
    print normal([array([1,0,0]), array([0,1,0]), array([1,1,0])])
    print normal([array([1,123,0123]), array([0345,1,034]), array([132,134,034])])
    print normal([array([1,0,0]), array([0,1,0]), array([1,123,1123])])