title: python-Interpolate polynomial where coefficients are matrices
body: <p>I have a polynomial of the form:</p>
  
  <pre><code>p(y) = A + By + Cy^2 ... + Dy^n
  </code></pre>
  
  <p>Here, each of the coefficients <code>A,B,..,D</code> are matrices (and therefore <code>p(y)</code> is also a matrix). Say I interpolate the polynomial at <code>n+1</code> points. I should now be able to solve this system. I'm trying to do this in Numpy. I have the following code right now:</p>
  
  <pre><code>a = np.vander([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2]) #polynomial degree is 12, a -&gt; (12x12)
  b = np.random.rand(12,60,60) #p(x) is a 60x60 matrix that I have evaluated at 12 points
  x = np.linalg.solve(a,b)
  </code></pre>
  
  <p>I get the following error:</p>
  
  <pre><code>ValueError: solve: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,m),(m,n)-&gt;(m,n) (size 60 is different from 12)
  </code></pre>
  
  <p>How can I solve this system in Numpy to get <code>x</code>? Is there a general mathematical trick to this?</p>