title: solving a sparse non linear system of equations using scipy.optimize.root
body: <p>I want to solve the following non-linear system of equations.</p>
  
  <p><a href="https://i.stack.imgur.com/EzY4U.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/EzY4U.png" alt="enter image description here"></a></p>
  
  <h3>Notes</h3>
  
  <ul>
  <li>the <code>dot</code> between <code>a_k</code> and <code>x</code> represents <code>dot product</code>.</li>
  <li>the <code>0</code> in the first equation represents <code>0 vector</code> and <code>0</code> in the second equation is <code>scaler 0</code></li>
  <li>all the matrices are sparse if that matters.</li>
  </ul>
  
  <h3>Known</h3>
  
  <ul>
  <li><code>K</code> is an <code>n x n</code> (positive definite) matrix</li>
  <li>each <code>A_k</code> is a known (symmetric) matrix</li>
  <li>each <code>a_k</code> is a known n x 1 vector</li>
  <li><code>N</code> is known (let's say N = 50). But I need a method where I can easily change N.</li>
  </ul>
  
  <h3>Unknown (trying to solve for)</h3>
  
  <ul>
  <li><code>x</code> is an <code>n x 1</code> a vector.</li>
  <li>each <code>alpha_k</code> for <code>1 &lt;= k &lt;= N</code>  a scaler</li>
  </ul>
  
  <h3>My thinking.</h3>
  
  <p>I am thinking of using <a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.html#multidimensional" rel="nofollow noreferrer">scipy root</a> to find x and each alpha_k. We essentially have <code>n</code> equations from each row of the first equation and another <code>N</code> equations from the constraint equations to solve for our <code>n + N</code> variables. Therefore we have the required number of equations to have a solution.</p>
  
  <p>I also have a reliable initial guess for <code>x</code> and the <code>alpha_k's</code>.</p>
  
  <h3>Toy example.</h3>
  
  <pre><code>n = 4
  N = 2
  K = np.matrix([[0.5, 0, 0, 0], [0, 1, 0, 0],[0,0,1,0], [0,0,0,0.5]])
  A_1 = np.matrix([[0.98,0,0.46,0.80],[0,0,0.56,0],[0.93,0.82,0,0.27],[0,0,0,0.23]])
  A_2 = np.matrix([[0.23, 0,0,0],[0.03,0.01,0,0],[0,0.32,0,0],[0.62,0,0,0.45]])
  a_1 = np.matrix(scipy.rand(4,1))
  a_2 = np.matrix(scipy.rand(4,1))
  </code></pre>
  
  <p>We are trying to solve for</p>
  
  <pre><code> x = [x1, x2, x3, x4] and alpha_1, alpha_2
  </code></pre>
  
  <h3>Questions:</h3>
  
  <ol>
  <li>I can actually brute force this toy problem and feed it to the solver. But how do I do I solve this toy problem in such a way that I can extend it easily to the case when I have let's say <code>n=50</code> and <code>N=50</code></li>
  <li>I will probably have to explicitly compute the Jacobian for larger matrices??.</li>
  </ol>
  
  <p>Can anyone give me any pointers?  </p>