How to find the 'centre' of vectors in 3-dimensional space.

Question

How do I find the 'centre' of a set of vectors in $ℝ^3$?

By centre I mean I want to find a vector such that maximum angle between it and any other vector in the set is kept to a minimum. Another way of looking at it would be to find the narrowest cone that contains all the vectors in the set.

I am currently working in MATLAB so methods that are simple to implement in code would be most welcome.

Answer 1

I assume that the vectors are contained within some half-space; otherwise any containing cone has to be greater than half-space which makes the geometry rather different.

First, normalize the vectors so that each of them has length $1$. They form a set $S=\{\vec u_j: j=1,\dots,N\}$ on the unit sphere. You want the smallest spherical cap that contains $S$. This is equivalent to finding a plane $P$ at maximal distance from $0$ that separates $S$ from $0$. Writing the plane as $P=\{\vec x: \vec x\cdot \vec w= 1\}$, we arrive at the convex minimization problem with linear constraints: $$\min \|\vec w\|^2 \quad \text{subject to } \ \vec u_j\cdot \vec w\ge 1,\quad j=1,\dots,N \tag1$$ The minimizing vector $\vec w$ is unique, and is the center vector that you wanted.

You can do all kinds of convex optimization in Matlab with CVX, and (1) is about the easiest convex optimization problem imaginable. Here it is in Sage:

vectors = [vector([1,2,3]), vector([1,0,2]), vector([3,2,4]), vector([5,2,-1]), vector([1,1,-1])] 
unit_vectors = [v/v.norm() for v in vectors]
constraints = [lambda x,vec=u: vector(x).dot_product(vec)-1 for u in unit_vectors]
target = lambda x: vector(x).norm()
minimize_constrained(target,constraints,[3,3,3])

Output: (1.38121931125,0.778255829438,0.427424333123), which looks reasonable.

Sage requires an initial point of optimization, which can be any vector that meets the constraints: you'll have it when you know a half-space containing your vectors.

---

I spent an inordinate amount of time pondering why the earlier version of this code didn't work, until I tried StackOverflow: Evaluating a list of python lambda functions only evaluates the last list element.

Answer 2

You could simply compute the mean of the vectors. Imagine a point cloud in $\mathbb{R}^3$ where each point is the end point of a vector starting at the origin. The mean vector is then a vector/point in the "middle" of this point cloud. In Matlab this is especially simple. Define a matrix $A$ with the vectors in its rows. If $N$ is the number of vectors, then $$\tt{m=sum(A)/N}$$ is the mean vector.