Programmatical Change of basis for coordinate vectors with different origin of coordinates (python/general maths)

Question

I have a Support Vector Machine that splits my data in two using a decision hyperplane (for visualisation purposes this is a sample dataset with three dimensions), like this: Now I want to perform a change of basis, such that the hyperplane lies flatly on the x/y plane, such that the distance from each sample point to the decision hyperplane is simply their z-coordinate.

For that, I know that I need to perform a change of basis. The hyperplane of the SVM is given by their coefficient (3d-vector) and intercept (scalar), using (as far as I understand it) the general form for mathematical planes: ax+by+cz=d, with a,b,c being the coordinates of the coefficient and d being the intercept. When plotted as 3d-Vector, the coefficient is a vector orthogonal to the plane (in the image it's the cyan line).

Now to the change of basis: If there was no intercept, I could just assume the vector that is the coefficient is one vector of my new basis, one other can be a random vector that is on the plane and the third one is simply cross product of both, resulting in three orthogonal vectors that can be the column vectors of the transformation-matrix. The z-function used in the code below comes from simple term rearrangement from the general form of planes: ax+by+cz=d <=> z=(d-ax-by)/c:

z_func = lambda interc, coef, x, y: (interc-coef[0]*x -coef[1]*y) / coef[2]
def generate_trafo_matrices(coefficient, z_func):
    normalize = lambda vec: vec/np.linalg.norm(vec)
    uvec2 = normalize(np.array([1, 0, z_func(1, 0)]))
    uvec3 = normalize(np.cross(uvec1, uvec2))
    back_trafo_matrix = np.array([uvec2, uvec3, coefficient]).T 
    #in other order such that its on the xy-plane instead of the yz-plane
    trafo_matrix = np.linalg.inv(back_trafo_matrix)
    return trafo_matrix, back_trafo_matrix

This transformation matrix would then be applied to all points, like this:

def _transform(self, points, inverse=False):
    trafo_mat = self.inverse_trafo_mat if inverse else self.trafo_mat              
    points = np.array([trafo_mat.dot(point) for point in points])        
    return points

Now if the intercept would be zero, that would work perfectly and the plane would be flat on the xy-axis. However as soon as I have an intercept != zero, the plane is not flat anymore:

I understand that that is the case because this is not a simple change of basis, because the coordinate origin of my other basis is not at (0,0,0) but at a different place (the hyperplane could be crossing the coefficient-vector at any point), but my attempts of adding the intercept to the transformation all didn't lead to the correct result:

def _transform(self, points, inverse=False):
    trafo_mat = self.inverse_trafo_mat if inverse else self.trafo_mat      
    intercept = self.intercept if inverse else -self.intercept
    ursprung_translate = trafo_mat.dot(np.array([0,0,0])+trafo_matrix[:,0]*intercept)
    points = np.array([point+trafo_matrix[:,0]*intercept for point in points])
    points = np.array([trafo_mat.dot(point) for point in points])        
    points = np.array([point-ursprung_translate for point in points])
    return points

is for example wrong. I am asking this on StackOverflow and not on the math StackExchange because I think I wouldn't be able to translate the respective math into code, I am glad I even got this far.

I have created a github gist with the code to do the transformation and create the plots at https://gist.github.com/cstenkamp/0fce4d662beb9e07f0878744c7214995, which can be launched using Binder under the link https://mybinder.org/v2/gist/jtpio/0fce4d662beb9e07f0878744c7214995/master?urlpath=lab%2Ftree%2Fchange_of_basis_with_translate.ipynb if somebody wants to play around with the code itself.

Any help is appreciated!

Answer 1

The problem here is that your plane is an affine space, not a vector space, so you can't use the usual transform matrix formula.

A coordinate system in affine space is given by an origin point and a basis (put together, they're called an affine frame). For example, if your origin is called O, the coordinates of the point M in the affine frame will be the cooordinates of the OM vector in the affine frame's basis.

As you can see, the "normal" R^3 space is a special case of affine space where the origin is (0,0,0).

Once we've determined those, we can use the frame change formulas in affine spaces: if we have two affine frames R = (O, b) and R' = (O', b'), the base change formula for a point M is: M(R') = base_change_matrix_from_b'_to_b * (M(R) - O'(R)) (with O'(R) the coordinates of O' in the coordinate system defined by R).

In our case, we're trying to go from the frame with an origin at (0,0,0) and the canonical basis, to a frame where the origin is the orthogonal projection of (0,0,0) on the plane and the basis is, for instance, the one described in your initial post.

Let's implement these steps:

To begin with, we'll define a Plane class to make our lifes a bit easier:

from dataclasses import dataclass
import numpy as np

@dataclass
class Plane:
    a: float
    b: float
    c: float
    d: float
    
    @property
    def normal(self):
        return np.array([self.a, self.b, self.c])
    
    def __contains__(self, point:np.array):
        return np.isclose(self.a*point[0] + self.b*point[1] + self.c*point[2] + self.d, 0)
    
    def project(self, point):
        x,y,z = point
        k = (self.a*x + self.b*y + self.c*z + self.d)/(self.a**2 + self.b**2 + self.c**2)
        return np.array([x - k*self.a, y-k*self.b, z-k*self.c])
   
    
    def z(self, x, y):
        return (- self.d - self.b*y - self.a*x)/self.c

We can then implement make_base_changer, which takes a Plane as an input, and return 2 lambda functions performing the forward and inverse transform (taking and returning a point). You should be able to test

def normalize(vec):
    return vec/np.linalg.norm(vec)
def make_base_changer(plane):
    uvec1 = plane.normal
    uvec2 = [0, -plane.d/plane.b, plane.d/plane.c]
    uvec3 = np.cross(uvec1, uvec2)
    transition_matrix = np.linalg.inv(np.array([uvec1, uvec2, uvec3]).T)
    
    origin = np.array([0,0,0])
    new_origin = plane.project(origin)
    forward  = lambda point: transition_matrix.dot(point - new_origin)
    backward = lambda point: np.linalg.inv(transition_matrix).dot(point) + new_origin
    return forward, backward