Solving a non-singular system of equations, hyperspace through n points in n-dim using Julia

Question

I have a simple problem with Julia: I would like to find the equation of an hyperplane passing through n points in dimension n. Is there any simple way to do that? I was solving a linear system of equations, however it can be non-singular and in this case, Julia returns an error. Is there any known way to solve parametric or non-singular systems of equations in Julia?

As an example, consider the 3d points [1 0 0], [0 0 1] and [1 0 1]. I would like to have y = 0 as a solution, via the vector [0 1 0 0] of coefficients.

In sage

a,b,c,d = var('a b c d')
f1 = a + d
f2 = a + c + d
f3 = c + d
solve([f1==0,f2==0,f3==0],a,b,c,d)

gives

[[a == 0, b == r1, c == 0, d == 0]]

Thank you very much for your help.

Answer 1

One way to do this is to find an eigenvector for the smallest eigenvalue of the covariance matrix of the points. This is a bit heavy computationally though.

We can represent a hyperplane by a unit vector n and a scalar d. The points of the hyperplane are then the x with n'*x+d = 0. For a point x n'*x+d is the distance of x from the plane. We can find the hyperplane by least squares:

If X[] are the k given points, we seek n, d to minimse

Q(n,d) = Sum{ i | sqr( X[i]'*n+d)}/k

Given n, the minimising d will be -Xbar.n, where Xbar is the mean of the X, so we may as well plug this d in and minimise

Q(n) = Sum{ i | sqr( x[i]'*n)}/k = n'*C*n 

where

x[i] = X[i] - Xbar
C = Sum{ i | x[i]*x[i]'}/k -- the covariance of the data

This Q will be minimised by an n that is an eigenvector for a minimal eigenvalue of C.

So the drill is: compute Xbar, and C. Diagonalise C (as U * D * U' say) and if a minimal entry of D is at (i,i), choose for n the i'th column of U), and finally compute d as -n'*Xbar.

Note that the only difficulty here is that there here may be many minimal values of D -- for example if all your points lay on a n-2 dimensional plane.

Answer 2

I'm not sure I completely understand the problem, but if we phrase above problem using a coefficient matrix and vector, Julia allows

julia> [1 0 0 1; 1 0 1 1; 0 0 1 1] \ [0; 0; 0]
4-element Array{Float64,1}:
 -0.0
  0.0
 -0.0
 -0.0

to get one particular solution and

julia> nullspace([1 0 0 1; 1 0 1 1; 0 0 1 1])
4×1 Array{Float64,2}:
  3.92523e-17
  1.0        
  0.0        
 -3.92523e-17

to get a base of the solution space. (Numerical issues here are a bit unfortunate, it should ideally be [0; 1; 0; 0], of course.) If we read that as [a; b; c; d] = [0; 0; 0; 0] + r*[0; 1; 0; 0], it's basically the solution given in the question.

If the matrix is rank-deficient, e.g. because the points are collinear, you will still get a particular solution from \, but nullspace will return a base with more that just one vector.