2D affine mapping using scipy

Question

I have two sets of 2D data A and B. I have to assume that there is a 2×2 matrix M and a 2x1 column vector c such that:

B = M*A + c + e

With e being the gaussian error (mean of 0 and variance of sigma^2).

My goal is to estimate the matrix M and vector c for which the sum of squared distances between B and the predictions is minimized.

So far, that's what I've done:

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

A = np.array([
    [0, -6],
    [0, 6],
    [-6, 0],
    [6, 0],
    [-3, -3],
    [-3, 3],
    [0,0],
    [2, -2],
    [2, 4],
    [-3, 0],
    [-5, -3],
    [2, 5]
])

B = np.array([
    [0.45070423, 0.52288732],
    [0.125, 0.5625],
    [0.1754386, 0.15789474],
    [0.200489, 0.55012225],
    [0.30769231, 0.32451923],
    [0.1375, 0.45],
    [0.19935691, 0.46302251],
    [0.32142857, 0.53571429],
    [0.12765957, 0.53191489],
    [0.15246637, 0.30269058],
    [0.2247191, 0.20224719],
    [0.14379085, 0.55882353],
])

# MAPPING : getting M and c

def error_function(params):
    # Unpack M and c from the flat parameter array
    M = params[:4].reshape(2, 2)
    c = params[4:]
    
    predicted_B = (A @ M) + c

    # Calculate the sum of squared distances
    return np.sum(np.sum((predicted_B - B) ** 2, axis=1))

# Perform the optimization
initial_guess = np.zeros(6)
result = minimize(error_function, initial_guess)
M = result.x[:4].reshape(2, 2)
c = result.x[4:]

# Getting predictions
predicted_B = (A @ M) + c

# Plotting
plt.scatter(B[:,0], B[:,1], color='red', s=5)
plt.scatter(predicted_B[:,0], predicted_B[:,1], color='blue', s=5)
plt.show()

Even if at the end of the optimization process I do get a low sum of squared distances (0.07), the result is not quite good if I plot the pairs of actual point/prediction on a graph.

Did I implement it correctly ? I am not 100% sure about the minimization of the sum of squared distances.

I guess also that the low sum of squared distances that I get is not that low since my targeted data are small values (between 0 and 0.7).

How could I improve it ? I have to stick strictly to the initial assumptions.

Answer 1

I'm assuming you want to solve (in the least squares sense) B = A @ M + c for M and c. This is a little different than what you wrote in the top post, but it's consistent with your code.

By inspection, your code appears to be correct for solving this problem, but that doesn't guarantee that it is correct, that minimize is able to solve the problem you pose, or that it cannot be improved.

To confirm that it was solving the problem correctly, my original approach was to use the analytical solution to B = A @ M within the objective function and to optimize for c only. With only two variables, it was easier to confirm manually that your solution was obtaining the correct result.

Upon further research, I see that if you form B2 = [B 1] (where 1 represents a column of 1s) and M2 = [[M], [x]], then B2 = A @ M2 is now of the form that can be solved analytically, and it is equivalent to your original problem.

So you can (very slightly) improve the accuracy of the solution by using the analytical solution, which can be computed with scipy.linalg.leastsq. After performing your imports and definition of A and B:

from scipy import linalg
A2 = np.concatenate((A, np.ones(len(A))[:, np.newaxis]), axis=1)
M2 = linalg.lstsq(A2, B)[0]
M2
# array([[ 0.00717165,  0.03502682],
#        [-0.02547716,  0.00720395],
#        [ 0.22712728,  0.45114461]])

You will find that these values match the ones produced by your optimization, but the solution is faster, accurate to more significant figures, and - most importantly - carries stronger guarantees of optimality.