Projection along affine hull

edit

Here is an example. It requires to load the Projection function:

sage: from sage.geometry.polyhedron.plot import Projection

Let's say we have the following polytope:

sage: P = Polyhedron(vertices=[[-2,0,3],[2,0,3],[0,2,1],[0,-2,1]])

and we are interested in the affine space spanned by the points

sage: affine_basis = [vector([1,1,2]),vector([1,-1,2]),vector([-1,-1,2])]

we first get a linear subspace and create the projection matrix:

sage: linear_subspace = [ap - affine_points[0] for ap in affine_basis[1:]]
sage: VS = VectorSpace(QQ,3)
sage: proj_matrix = matrix([sum([v.inner_product(lv)*lv/(lv.inner_product(lv)) for lv in linear_subspace]) for v in VS.basis()]).transpose()

Then, we create a function that will emulate the affine projection (notice the addition of the first element of the affine basis).

sage: def my_proj(x): return proj_matrix*x + affine_basis[0]

Then, we project the polytope and use the transformed coordinates of the vertices to create a new polytope.

sage: proj_p = Projection(P,proj=my_proj)
sage: projected_P = Polyhedron(vertices=proj_p.transformed_coords)

Now, we can see the difference between

sage: P
A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices
sage: P.vertices()
(A vertex at (-2, 0, 3),
 A vertex at (0, -2, 1),
 A vertex at (0, 2, 1),
 A vertex at (2, 0, 3))

and the resulting polytope:

sage: projected_P
A 2-dimensional polyhedron in QQ^3 defined as the convex hull of 4 vertices
sage: projected_P.vertices()
(A vertex at (0, -2, 2),
 A vertex at (0, 0, 2),
 A vertex at (2, 2, 2),
 A vertex at (2, 4, 2))

In sage 9.1+ matrices act on polyhedra as linear transformation. This makes it easy.

Let me explain how to project the 3-dimensional permutahedron P to the affine span of [(3, 2, 4, 1), (4, 1, 3, 2), (1, 2, 3, 4)]. By translation, we can reduce the problem to a linear one. We can instead linearly project P1 = P - vector((3, 2, 4, 1)) to the linear span of [((1, -1, -1, 1), (-2, 0, -1, 3))].

Obtaining P1:

sage: P = polytopes.permutahedron(4)  # dim 3 and ambient dim 4
sage: P1 = P - vector((3, 2, 4, 1))

To do this we use the right kernel:

sage: M = Matrix(((1, -1, -1, 1), (-2, 0, -1, 3)))
sage: B = M.stack(M.right_kernel_matrix()).transpose()
sage: B
[ 1 -2  1  0]
[-1  0  1  2]
[-1 -1  1 -3]
[ 1  3  1 -1]

The first two rows are spanning our linear space we want to map to. The third and fourth row are a base extension such that they are orthogonal to the first and second row.

Now we want to linearly transform P such that the vectors corresponding to the first and second row get mapped to themselves and the others to zero. That's linear algebra:

sage: P2 = B*Matrix([[1,0,0,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]])*B.inverse()*P1

Finally, we add the vector again:

sage: P2 + vector((3, 2, 4, 1))
A 2-dimensional polyhedron in QQ^4 defined as the convex hull of 10 vertices