Least-squares

For fun this week, I’ll use Jupyter notebooks inside Sage, which you create like this:

To turn on pretty output (typesetting), type %display typeset.

As far as I know, Sage does not have a built-in method to find a “least-squares solution” to a system of linear equations. The description of a least squares solution to Ax=b as a solution to A^TAx=A^Tb is easy to work with in Sage. For example, to find a least-squares solution to the system

x₁ + x₂ = 3
2x₁ + 3x₂ = 1
x₁ – x₂ = 5

we can type:
A = matrix([[1, 1],[2,3],[1,-1]]); A
B = A.transpose(); B
v = vector([3,1,5]); v
(B*A).solve_right(B*v)

or, more succinctly,
A = matrix([[1, 1],[2,3],[1,-1]])
(A.transpose()*A).solve_right(A.transpose()*vector([3,1,5]))

We could also write a function that finds least-squares solutions for us:

def least_squares(A,v):
     return (A.transpose()*A).solve_right(A.transpose()*v)

(Again, indentation is important: this is Python.) Then we can call it with the matrix and vector of your choice:
B = matrix([[1,2],[3,4],[5,6]])
v = vector([7,5,9])
B.solve_right(v)
least_squares(B,v)

(There is also built-in code in Sage and numpy for fitting curves to data.)