Geometric interpretation of Least Squares solution in an inconsistent system

311 Views Asked by At

Let's say we have an inconsistent set of three linear equations in two variables. We can interpret these as three lines forming a triangle in the plane:

Let us write this in matrix form as

$$A \begin{bmatrix}x \\ y\end{bmatrix} = b$$

If we now "solve" this system by determining

$$ \begin{bmatrix} \hat x \\ \hat y \end{bmatrix} = \arg\min \left\Vert A \begin{bmatrix}x \\ y\end{bmatrix} - b \right\Vert_2^2$$

is there a geometric interpretation of the point $(\hat x, \hat y)$ with respect to the triangle?

If so, is there a way to extend this notion to higher dimensions and/or more equations?


EDIT: I noticed that my quesiton was ill posed: You can scale any of the three equations with some nonzero factor and get a different minimizer so in that sense my original question was ill posed. Let us now assume that all the equations are normalized, that is, that each row vector of the matrix $A$ has norm $1$.

1

There are 1 best solutions below

3
On BEST ANSWER

If the rows of $A$ are normalized to 1, then the expression $$ \left( A_i \left[ x \atop y\right] - b_i\right)^2,$$ where $A_i$ is the $i$th row of $A$, denotes the distance of the point $(x,y)$ to the $i$th face of the polygon. Then, $$\left\| A \left[ x \atop y\right] -b \right\|^2$$ is the sum of the squared distances of the point $(x,y)$ to all faces of the polygon and you are looking for the point that minimizes this sum of distances.

In two dimensions, and with three lines, that would be the symmedian, or Lemoine, point. symmedian point

Code to generate the figure


import numpy as np
import matplotlib.pyplot as plt

A = np.asarray([
  [2.0, 1.0],
  [-0.5, 1.0],
  [-3.0, 1.0]
])

b = np.asarray([-1.0,1.0,-2.0])

for i in range(np.size(A,0)):
  print(A[i,:])
  print(np.linalg.norm(A[i,:]))
  A[i,:] = A[i,:]/np.linalg.norm(A[i,:])
  print(A[i,:])

print(A)


def line(i, x):
  return (-b[i] + x*A[i,0])/A[i,1]

if __name__ == "__main__":

  x0 = -1.5
  x1 = 3

  plt.plot([x0, x1], [line(0,x0), line(0,x1)])
  plt.plot([x0, x1], [line(1,x0), line(1,x1)])
  plt.plot([x0, x1], [line(2,x0), line(2,x1)])

  p_hat = np.dot(np.linalg.inv(np.dot(A.T,A)), np.dot(A.T, b))

  plt.scatter(p_hat[0], p_hat[1])
  plt.xlim(-1.5,5)
  plt.ylim(-3,4)
  plt.gca().set_aspect('equal')
  plt.show()