title: Modified Least Squares Fit in python with unique errors, solving for a matrix
body: <p>I have an over-determined system I want to minimize, very similar to (but not quite the same as) linear least squares. It is of this form:</p>
  <p>AX = Y</p>
  <p>where A will need to be 2x3, because X is 3xN and Y is 2xN.</p>
  <p>Here is the issue: the unknown here is the transformation matrix A, not X in the usual formulation of least squares. (also, I would mention that X and Y are matrices in this case, where they are usually vectors in classic least squares)</p>
  <p>Further complicating things, I want this to be a weighted least squares, where for each 2x1 vector in Y, (Y is 2xN, a stack of these vectors), I have a 2x1 weight vector that element-wise weights the vectors in Y (like np.divide(Y[i], np.square(weights[i]))).</p>
  <p>I first tried to solve things without the weights, and I got about this far:</p>
  <p>AX = Y</p>
  <p>X.T * A.T = Y.T</p>
  <p>X * X.T * A.T = X * Y.T</p>
  <p>A.T = (X * X.T)^-1 * X * Y.T</p>
  <p>A = Y * X.T * (X * X.T)^-1</p>
  <p>This resulted in this Python function:</p>
  <pre><code>def mat_solve(X,Y):
  return Y@X.T@np.linalg.inv(X@X.T)
  </code></pre>
  <p>But I can't obviously see that this truly minimizes the cost function of the problem, so I am not confident its truly right, even without the weights.</p>
  <p>I also tried writing the cost function out as a sum, and taking a derivative and setting it equal to 0, which resulted in this Python function:</p>
  <pre><code>def sum_solve(X,Y):
  result1 = 0
  result2 = 0
  for i in range(len(X)):
      result1 += Y[:,i:i+1]@X.T[i:i+1,:]
      result2 += X[:,i:i+1]@X.T[i:i+1,:]
  return result1@np.linalg.inv(result2)
  </code></pre>
  <p>These two functions return a very similar matrix, but its not quite the same. It should be the same if they both are truly solutions to the problem. I would mention that the space that is being minimized over is NOT nonlinear, and that I am very confident there is only one minimum.</p>
  <p>I thought I would go a more brute force route, and try using an autodiff package called JAX which has its own minimize() function very similar to scipy (but in a gradient-descent sort of way), and just write the cost function out explicitly, even with the weights in there. The weights were just turned into a stack of identity matrices and then multiplied by their Y[i]s, as you can see (I called the cost function myLSQ, and then just called minimize() on it, so it gradient descents with respect to A):</p>
  <pre><code>import jax 
  import jax.numpy as jnp
  from jax.scipy.optimize import minimize
  
  uncertainty_new = [np.identity(2)/(x**2) for x in uncertainty_list]
  freq_list_new = [np.append(x,[1]) for x in freq_list]
  
  def myLSQ(A,X,Y,W):
      Anew = A.reshape((2,3))
      temp = []
      result = 0
      for i in range(len(X)):
          temp = jnp.subtract(jnp.dot(Anew,X[i]), Y[i])
          result += jnp.dot(jnp.dot(temp,W[i]), temp.T)
      return result
  
  
  
  A_guess = np.array([[1.,1.],[1.,1.],[1.,1.]])
  %prun result = minimize(myLSQ, A_guess.flatten(),method = &quot;BFGS&quot;,args = (freq_list_new, beam_peak_list, uncertainty_new,))
  </code></pre>
  <p>But this thing runs for more than 20 min, and then returns an absolutely horrendous matrix that gets no where close to minimizing the problem.</p>
  <p>I apologize for the enormous essay, but here are my questions.</p>
  <ol>
  <li><p>Why do my two (unweighted) methods give slightly different results, and how do I add the weights in to them?</p>
  </li>
  <li><p>Is there a better way to do this? If so, what is it?</p>
  </li>
  <li><p>I want the fit to be weighted, but I don't want to have to give the A matrix weights when I use it later (once its correct). I want to be able to just do AX = Y, without weights on Y. The weights are just to modify the fit. How do I do this?</p>
  </li>
  </ol>