code: def MatrixSolveLs(a, rhs, l2_reg):
      
      r = np.empty(rhs.shape).astype(a.dtype)
      for coord in np.ndindex(a.shape[:-2]):
          pos = coord + (Ellipsis,)
          r[pos] = np.linalg.lstsq(a[pos], rhs[pos])[0]
  
      return r,
docs: Matrix least-squares solve op.