<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Sparse Linear Equations — CVXOPT User's Guide</title>
<link rel="stylesheet" href="_static/cvxopt.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: '#',
VERSION: '1.1.1',
COLLAPSE_MODINDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: false
};
</script>
<script type="text/javascript" src="_static/jquery.js"></script>
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="copyright" title="Copyright" href="copyright.html" />
<link rel="top" title="CVXOPT User's Guide" href="index.html" />
<link rel="next" title="Cone Programming" href="coneprog.html" />
<link rel="prev" title="Discrete Transforms" href="fftw.html" />
</head>
<body>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="coneprog.html" title="Cone Programming"
accesskey="N">next</a></li>
<li class="right" >
<a href="fftw.html" title="Discrete Transforms"
accesskey="P">previous</a> |</li>
<li><a href="http://abel.ee.ucla.edu/cvxopt">CVXOPT home</a> |</li>
<li><a href="index.html">user's guide</a> </li>
</ul>
</div>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body">
<div class="section" id="sparse-linear-equations">
<span id="c-spsolvers"></span><h1>Sparse Linear Equations<a class="headerlink" href="#sparse-linear-equations" title="Permalink to this headline">¶</a></h1>
<p>In this section we describe routines for solving sparse sets of linear
equations.</p>
<p>A real symmetric or complex Hermitian sparse matrix is stored as an
<a title="cvxopt.spmatrix" class="reference external" href="matrices.html#cvxopt.spmatrix"><tt class="xref docutils literal"><span class="pre">spmatrix</span></tt></a> object <tt class="docutils literal"><span class="pre">X</span></tt> of size
(<img class="math" src="_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/>, <img class="math" src="_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/>) and an
additional character argument <tt class="docutils literal"><span class="pre">uplo</span></tt> with possible values <tt class="xref docutils literal"><span class="pre">'L'</span></tt>
and <tt class="xref docutils literal"><span class="pre">'U'</span></tt>. If <tt class="docutils literal"><span class="pre">uplo</span></tt> is <tt class="xref docutils literal"><span class="pre">'L'</span></tt>, the lower triangular part
of <tt class="docutils literal"><span class="pre">X</span></tt> contains the lower triangular part of the symmetric or Hermitian
matrix, and the upper triangular matrix of <tt class="docutils literal"><span class="pre">X</span></tt> is ignored. If <tt class="docutils literal"><span class="pre">uplo</span></tt>
is <tt class="xref docutils literal"><span class="pre">'U'</span></tt>, the upper triangular part of <tt class="docutils literal"><span class="pre">X</span></tt> contains the upper
triangular part of the matrix, and the lower triangular matrix of <tt class="docutils literal"><span class="pre">X</span></tt> is
ignored.</p>
<p>A general sparse square matrix of order <img class="math" src="_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/> is represented by an
<tt class="xref docutils literal"><span class="pre">spmatrix</span></tt> object of size (<img class="math" src="_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/>, <img class="math" src="_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/>).</p>
<p>Dense matrices, which appear as righthand sides of equations, are
stored using the same conventions as in the BLAS and LAPACK modules.</p>
<div class="section" id="matrix-orderings">
<span id="s-orderings"></span><h2>Matrix Orderings<a class="headerlink" href="#matrix-orderings" title="Permalink to this headline">¶</a></h2>
<p>CVXOPT includes an interface to the AMD library for computing approximate
minimum degree orderings of sparse matrices.</p>
<div class="admonition-see-also admonition seealso">
<p class="first admonition-title">See also</p>
<ul class="last simple">
<li><a class="reference external" href="http://www.cise.ufl.edu/research/sparse/amd">AMD code, documentation, copyright, and license</a></li>
<li>P. R. Amestoy, T. A. Davis, I. S. Duff, Algorithm 837: AMD, An
Approximate Minimum Degree Ordering Algorithm, ACM Transactions on
Mathematical Software, 30(3), 381-388, 2004.</li>
</ul>
</div>
<dl class="function">
<dt id="cvxopt.amd.order">
<tt class="descclassname">cvxopt.amd.</tt><tt class="descname">order</tt><big>(</big><em>A</em><span class="optional">[</span>, <em>uplo = 'L'</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.amd.order" title="Permalink to this definition">¶</a></dt>
<dd>Computes the approximate mimimum degree ordering of a symmetric sparse
matrix <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/>. The ordering is returned as an integer dense matrix
with length equal to the order of <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/>. Its entries specify a
permutation that reduces fill-in during the Cholesky factorization.
More precisely, if <tt class="docutils literal"><span class="pre">p</span> <span class="pre">=</span> <span class="pre">order(A)</span></tt> , then <tt class="docutils literal"><span class="pre">A[p,</span> <span class="pre">p]</span></tt> has
sparser Cholesky factors than <tt class="docutils literal"><span class="pre">A</span></tt>.</dd></dl>
<p>As an example we consider the matrix</p>
<div class="math">
<p><img src="_images/math/d6eaef0ba36a702d7918e25eeb5b544e3f24f035.png" alt="\left[ \begin{array}{rrrr}
10 & 0 & 3 & 0 \\
0 & 5 & 0 & -2 \\
3 & 0 & 5 & 0 \\
0 & -2 & 0 & 2
\end{array}\right]." /></p>
</div><div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">amd</span>
<span class="gp">>>> </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">([</span><span class="mi">10</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="o">-</span><span class="mi">2</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">])</span>
<span class="gp">>>> </span><span class="n">P</span> <span class="o">=</span> <span class="n">amd</span><span class="o">.</span><span class="n">order</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">P</span>
<span class="go">[ 1]</span>
<span class="go">[ 0]</span>
<span class="go">[ 2]</span>
<span class="go">[ 3]</span>
</pre></div>
</div>
</div>
<div class="section" id="general-linear-equations">
<span id="s-umfpack"></span><h2>General Linear Equations<a class="headerlink" href="#general-linear-equations" title="Permalink to this headline">¶</a></h2>
<p>The module <tt class="xref docutils literal"><span class="pre">cvxopt.umfpack</span></tt> includes four functions for solving
sparse non-symmetric sets of linear equations. They call routines from
the UMFPACK library, with all control options set to the default values
described in the UMFPACK user guide.</p>
<div class="admonition-see-also admonition seealso">
<p class="first admonition-title">See also</p>
<ul class="last simple">
<li><a class="reference external" href="http://www.cise.ufl.edu/research/sparse/umfpack">UMFPACK code, documentation, copyright, and license</a></li>
<li>T. A. Davis, Algorithm 832: UMFPACK – an unsymmetric-pattern
multifrontal method with a column pre-ordering strategy, ACM
Transactions on Mathematical Software, 30(2), 196-199, 2004.</li>
</ul>
</div>
<dl class="function">
<dt id="cvxopt.umfpack.linsolve">
<tt class="descclassname">cvxopt.umfpack.</tt><tt class="descname">linsolve</tt><big>(</big><em>A</em>, <em>B</em><span class="optional">[</span>, <em>trans = 'N'</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.umfpack.linsolve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves a sparse set of linear equations</p>
<div class="math">
<p><img src="_images/math/ff3e601d603ae2c0fef98fe8b78d53fb099c9ee9.png" alt="AX & = B \quad (\mathrm{trans} = \mathrm{'N'}), \\
A^TX & = B \quad (\mathrm{trans} = \mathrm{'T'}), \\
A^HX & = B \quad (\mathrm{trans} = \mathrm{'C'})," /></p>
</div><p>where <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/> is a sparse matrix and <img class="math" src="_images/math/ff5fb3d775862e2123b007eb4373ff6cc1a34d4e.png" alt="B"/> is a dense matrix.
The arguments <tt class="docutils literal"><span class="pre">A</span></tt> and <tt class="docutils literal"><span class="pre">B</span></tt> must have the same type
(<tt class="xref docutils literal"><span class="pre">'d'</span></tt> or <tt class="xref docutils literal"><span class="pre">'z'</span></tt>) as <tt class="docutils literal"><span class="pre">A</span></tt>. On exit <tt class="docutils literal"><span class="pre">B</span></tt> contains
the solution. Raises an <tt class="xref docutils literal"><span class="pre">ArithmeticError</span></tt> if the coefficient
matrix is singular.</p>
</dd></dl>
<p>In the following example we solve an equation with coefficient matrix</p>
<div class="math" id="equation-e-sp-Adef">
<p><span class="eqno">(1)</span><img src="_images/math/6fea0b5ad40e6222f621a260aeb1fc3d0b2d1e59.png" alt="A = \left[\begin{array}{rrrrr}
2 & 3 & 0 & 0 & 0 \\
3 & 0 & 4 & 0 & 6 \\
0 &-1 &-3 & 2 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 4 & 2 & 0 & 1
\end{array}\right]." /></p>
</div><div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">umfpack</span>
<span class="gp">>>> </span><span class="n">V</span> <span class="o">=</span> <span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">I</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">J</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="n">V</span><span class="p">,</span><span class="n">I</span><span class="p">,</span><span class="n">J</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">B</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">5</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
<span class="gp">>>> </span><span class="n">umfpack</span><span class="o">.</span><span class="n">linsolve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="n">B</span><span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">B</span>
<span class="go">[ 5.79e-01]</span>
<span class="go">[-5.26e-02]</span>
<span class="go">[ 1.00e+00]</span>
<span class="go">[ 1.97e+00]</span>
<span class="go">[-7.89e-01]</span>
</pre></div>
</div>
<p>The function <a title="cvxopt.umfpack.linsolve" class="reference internal" href="#cvxopt.umfpack.linsolve"><tt class="xref docutils literal"><span class="pre">linsolve</span></tt></a> is
equivalent to the following three functions called in sequence.</p>
<dl class="function">
<dt id="cvxopt.umfpack.symbolic">
<tt class="descclassname">cvxopt.umfpack.</tt><tt class="descname">symbolic</tt><big>(</big><em>A</em><big>)</big><a class="headerlink" href="#cvxopt.umfpack.symbolic" title="Permalink to this definition">¶</a></dt>
<dd>Reorders the columns of <tt class="docutils literal"><span class="pre">A</span></tt> to reduce fill-in and performs a symbolic
LU factorization. <tt class="docutils literal"><span class="pre">A</span></tt> is a sparse, possibly rectangular, matrix.
Returns the symbolic factorization as an opaque C object that can be
passed on to <a title="cvxopt.umfpack.numeric" class="reference internal" href="#cvxopt.umfpack.numeric"><tt class="xref docutils literal"><span class="pre">numeric</span></tt></a>.</dd></dl>
<dl class="function">
<dt id="cvxopt.umfpack.numeric">
<tt class="descclassname">cvxopt.umfpack.</tt><tt class="descname">numeric</tt><big>(</big><em>A</em>, <em>F</em><big>)</big><a class="headerlink" href="#cvxopt.umfpack.numeric" title="Permalink to this definition">¶</a></dt>
<dd>Performs a numeric LU factorization of a sparse, possibly rectangular,
matrix <tt class="docutils literal"><span class="pre">A</span></tt>. The argument <tt class="docutils literal"><span class="pre">F</span></tt> is the symbolic factorization
computed by <a title="cvxopt.umfpack.symbolic" class="reference internal" href="#cvxopt.umfpack.symbolic"><tt class="xref docutils literal"><span class="pre">symbolic</span></tt></a>
applied to the matrix <tt class="docutils literal"><span class="pre">A</span></tt>,
or another sparse matrix with the same sparsity pattern, dimensions,
and type. The numeric factorization is returned as an opaque C object
that that can be passed on to
<a title="cvxopt.umfpack.solve" class="reference internal" href="#cvxopt.umfpack.solve"><tt class="xref docutils literal"><span class="pre">solve</span></tt></a>. Raises an
<tt class="xref docutils literal"><span class="pre">ArithmeticError</span></tt> if the matrix is singular.</dd></dl>
<dl class="function">
<dt id="cvxopt.umfpack.solve">
<tt class="descclassname">cvxopt.umfpack.</tt><tt class="descname">solve</tt><big>(</big><em>A</em>, <em>F</em>, <em>B</em><span class="optional">[</span>, <em>trans = 'N'</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.umfpack.solve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves a set of linear equations</p>
<div class="math">
<p><img src="_images/math/ff3e601d603ae2c0fef98fe8b78d53fb099c9ee9.png" alt="AX & = B \quad (\mathrm{trans} = \mathrm{'N'}), \\
A^TX & = B \quad (\mathrm{trans} = \mathrm{'T'}), \\
A^HX & = B \quad (\mathrm{trans} = \mathrm{'C'})," /></p>
</div><p>where <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/> is a sparse matrix and <img class="math" src="_images/math/ff5fb3d775862e2123b007eb4373ff6cc1a34d4e.png" alt="B"/> is a dense matrix.
The arguments <tt class="docutils literal"><span class="pre">A</span></tt> and <tt class="docutils literal"><span class="pre">B</span></tt> must have the same type. The argument
<tt class="docutils literal"><span class="pre">F</span></tt> is a numeric factorization computed
by <a title="cvxopt.umfpack.numeric" class="reference internal" href="#cvxopt.umfpack.numeric"><tt class="xref docutils literal"><span class="pre">numeric</span></tt></a>.
On exit <tt class="docutils literal"><span class="pre">B</span></tt> is overwritten by the
solution.</p>
</dd></dl>
<p>These separate functions are useful for solving several sets of linear
equations with the same coefficient matrix and different righthand sides,
or with coefficient matrices that share the same sparsity pattern.
The symbolic factorization depends only on the sparsity pattern of
the matrix, and not on the numerical values of the nonzero coefficients.
The numerical factorization on the other hand depends on the sparsity
pattern of the matrix and on its the numerical values.</p>
<p>As an example, suppose <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/> is the matrix <a href="#equation-e-sp-Adef">(1)</a> and</p>
<div class="math">
<p><img src="_images/math/df54abf5204a3ea7bb2efce847495127e99154e4.png" alt="B = \left[\begin{array}{rrrrr}
4 & 3 & 0 & 0 & 0 \\
3 & 0 & 4 & 0 & 6 \\
0 &-1 &-3 & 2 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 4 & 2 & 0 & 2
\end{array}\right]," /></p>
</div><p>which differs from <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/> in its first and last entries. The following
code computes</p>
<div class="math">
<p><img src="_images/math/4f017986e2e23bd929c17c8463200dea3fa870df.png" alt="\newcommand{\ones}{\mathbf 1}
x = A^{-T}B^{-1}A^{-1}\ones." /></p>
</div><div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">umfpack</span>
<span class="gp">>>> </span><span class="n">VA</span> <span class="o">=</span> <span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">VB</span> <span class="o">=</span> <span class="p">[</span><span class="mi">4</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">2</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">I</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">J</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="n">VA</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">B</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="n">VB</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">x</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">5</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
<span class="gp">>>> </span><span class="n">Fs</span> <span class="o">=</span> <span class="n">umfpack</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">FA</span> <span class="o">=</span> <span class="n">umfpack</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">Fs</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">FB</span> <span class="o">=</span> <span class="n">umfpack</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">B</span><span class="p">,</span> <span class="n">Fs</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">umfpack</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">FA</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">umfpack</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">B</span><span class="p">,</span> <span class="n">FB</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">umfpack</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">FA</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">trans</span><span class="o">=</span><span class="s">'T'</span><span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">x</span>
<span class="go">[ 5.81e-01]</span>
<span class="go">[-2.37e-01]</span>
<span class="go">[ 1.63e+00]</span>
<span class="go">[ 8.07e+00]</span>
<span class="go">[-1.31e-01]</span>
</pre></div>
</div>
</div>
<div class="section" id="positive-definite-linear-equations">
<span id="s-cholmod"></span><h2>Positive Definite Linear Equations<a class="headerlink" href="#positive-definite-linear-equations" title="Permalink to this headline">¶</a></h2>
<p><tt class="xref docutils literal"><span class="pre">cvxopt.cholmod</span></tt> is an interface to the Cholesky factorization routines
of the CHOLMOD package. It includes functions for Cholesky factorization
of sparse positive definite matrices, and for solving sparse sets of linear
equations with positive definite matrices.
The routines can also be used for computing
<span class="raw-html">LDL<sup><small>T</small></sup></span>
(or
<span class="raw-html">LDL<sup><small>H</small></sup></span>
factorizations
of symmetric indefinite matrices (with <img class="math" src="_images/math/859ccf4cd60c7bc6b8fa1afc9a42dc811a826d6f.png" alt="L"/> unit lower-triangular and
<img class="math" src="_images/math/9ffb448918db29f2a72f8f87f421b3b3cad18f95.png" alt="D"/> diagonal and nonsingular) if such a factorization exists.</p>
<div class="admonition-see-also admonition seealso">
<p class="first admonition-title">See also</p>
<p class="last"><a class="reference external" href="http://www.cise.ufl.edu/research/sparse/cholmod">CHOLMOD code, documentation, copyright, and license</a></p>
</div>
<dl class="function">
<dt id="cvxopt.cholmod.linsolve">
<tt class="descclassname">cvxopt.cholmod.</tt><tt class="descname">linsolve</tt><big>(</big><em>A</em>, <em>B</em><span class="optional">[</span>, <em>p = None</em>, <em>uplo = 'L'</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.cholmod.linsolve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves</p>
<div class="math">
<p><img src="_images/math/57a233e4bced8a9e14812ddb76241c0e78d32d86.png" alt="AX = B" /></p>
</div><p>with <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/> sparse and real symmetric or complex Hermitian.</p>
<p><tt class="docutils literal"><span class="pre">B</span></tt> is a dense matrix of the same type as <tt class="docutils literal"><span class="pre">A</span></tt>. On exit it
is overwritten with the solution. The argument <tt class="docutils literal"><span class="pre">p</span></tt> is an integer
matrix with length equal to the order of <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/>, and specifies an
optional reordering. If <tt class="docutils literal"><span class="pre">p</span></tt> is not specified, CHOLMOD uses a
reordering from the AMD library.</p>
<p>Raises an <tt class="xref docutils literal"><span class="pre">ArithmeticError</span></tt> if the factorization does not exist.</p>
</dd></dl>
<p>As an example, we solve</p>
<div class="math" id="equation-e-A-pd">
<p><span class="eqno">(2)</span><img src="_images/math/5589a714a734f3834cc4aead152691d7517b01d0.png" alt="\left[ \begin{array}{rrrr}
10 & 0 & 3 & 0 \\
0 & 5 & 0 & -2 \\
3 & 0 & 5 & 0 \\
0 & -2 & 0 & 2
\end{array}\right] X =
\left[ \begin{array}{cc}
0 & 4 \\ 1 & 5 \\ 2 & 6 \\ 3 & 7
\end{array} \right]." /></p>
</div><div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">cholmod</span>
<span class="gp">>>> </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">([</span><span class="mi">10</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mi">2</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">])</span>
<span class="gp">>>> </span><span class="n">X</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">8</span><span class="p">),</span> <span class="p">(</span><span class="mi">4</span><span class="p">,</span><span class="mi">2</span><span class="p">),</span> <span class="s">'d'</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">cholmod</span><span class="o">.</span><span class="n">linsolve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="n">X</span><span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">X</span>
<span class="go">[-1.46e-01 4.88e-02]</span>
<span class="go">[ 1.33e+00 4.00e+00]</span>
<span class="go">[ 4.88e-01 1.17e+00]</span>
<span class="go">[ 2.83e+00 7.50e+00]</span>
</pre></div>
</div>
<dl class="function">
<dt id="cvxopt.cholmod.splinsolve">
<tt class="descclassname">cvxopt.cholmod.</tt><tt class="descname">splinsolve</tt><big>(</big><em>A</em>, <em>B</em><span class="optional">[</span>, <em>p = None</em>, <em>uplo = 'L'</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.cholmod.splinsolve" title="Permalink to this definition">¶</a></dt>
<dd>Similar to <a title="cvxopt.cholmod.linsolve" class="reference internal" href="#cvxopt.cholmod.linsolve"><tt class="xref docutils literal"><span class="pre">linsolve</span></tt></a> except that
<tt class="docutils literal"><span class="pre">B</span></tt> is an <a title="cvxopt.spmatrix" class="reference external" href="matrices.html#cvxopt.spmatrix"><tt class="xref docutils literal"><span class="pre">spmatrix</span></tt></a> and
that the solution is returned as an output argument (as a new
<tt class="xref docutils literal"><span class="pre">spmatrix</span></tt>). <tt class="docutils literal"><span class="pre">B</span></tt> is not modified.</dd></dl>
<p>The following code computes the inverse of the coefficient matrix
in <a href="#equation-e-A-pd">(2)</a> as a sparse matrix.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">X</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">splinsolve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">spmatrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span><span class="nb">range</span><span class="p">(</span><span class="mi">4</span><span class="p">),</span><span class="nb">range</span><span class="p">(</span><span class="mi">4</span><span class="p">)))</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">X</span>
<span class="go">[ 1.22e-01 0 -7.32e-02 0 ]</span>
<span class="go">[ 0 3.33e-01 0 3.33e-01]</span>
<span class="go">[-7.32e-02 0 2.44e-01 0 ]</span>
<span class="go">[ 0 3.33e-01 0 8.33e-01]</span>
</pre></div>
</div>
<p>The functions <a title="cvxopt.cholmod.linsolve" class="reference internal" href="#cvxopt.cholmod.linsolve"><tt class="xref docutils literal"><span class="pre">linsolve</span></tt></a> and
<a title="cvxopt.cholmod.splinsolve" class="reference internal" href="#cvxopt.cholmod.splinsolve"><tt class="xref docutils literal"><span class="pre">splinsolve</span></tt></a> are equivalent to
<a title="cvxopt.cholmod.symbolic" class="reference internal" href="#cvxopt.cholmod.symbolic"><tt class="xref docutils literal"><span class="pre">symbolic</span></tt></a> and
<a title="cvxopt.cholmod.numeric" class="reference internal" href="#cvxopt.cholmod.numeric"><tt class="xref docutils literal"><span class="pre">numeric</span></tt></a> called in sequence, followed by
<a title="cvxopt.cholmod.solve" class="reference internal" href="#cvxopt.cholmod.solve"><tt class="xref docutils literal"><span class="pre">solve</span></tt></a>, respectively,
<a title="cvxopt.cholmod.spsolve" class="reference internal" href="#cvxopt.cholmod.spsolve"><tt class="xref docutils literal"><span class="pre">spsolve</span></tt></a>.</p>
<dl class="function">
<dt id="cvxopt.cholmod.symbolic">
<tt class="descclassname">cvxopt.cholmod.</tt><tt class="descname">symbolic</tt><big>(</big><em>A</em><span class="optional">[</span>, <em>p = None</em>, <em>uplo = 'L'</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.cholmod.symbolic" title="Permalink to this definition">¶</a></dt>
<dd><p>Performs a symbolic analysis of a sparse real symmetric or
complex Hermitian matrix <img class="math" src="_images/math/019e9892786e493964e145e7c5cf7b700314e53b.png" alt="A"/> for one of the two factorizations:</p>
<div class="math" id="equation-e-chol-ll">
<p><span class="eqno">(3)</span><img src="_images/math/7828b34252e1af5b450602189daf8af2e6217c2a.png" alt="PAP^T = LL^T, \qquad PAP^T = LL^H," /></p>
</div><p>and</p>
<div class="math" id="equation-e-chol-ldl">
<p><span class="eqno">(4)</span><img src="_images/math/1c839f56dd2ed41d4e9d5b2f84296b1c561dfd8b.png" alt="PAP^T = LDL^T, \qquad PAP^T = LDL^H," /></p>
</div><p>where <img class="math" src="_images/math/4b4cade9ca8a2c8311fafcf040bc5b15ca507f52.png" alt="P"/> is a permutation matrix, <img class="math" src="_images/math/859ccf4cd60c7bc6b8fa1afc9a42dc811a826d6f.png" alt="L"/> is lower triangular
(unit lower triangular in the second factorization), and <img class="math" src="_images/math/9ffb448918db29f2a72f8f87f421b3b3cad18f95.png" alt="D"/> is
nonsingular diagonal. The type of factorization depends on the value
of <tt class="xref docutils literal"><span class="pre">options['supernodal']</span></tt> (see below).</p>
<p>If <tt class="docutils literal"><span class="pre">uplo</span></tt> is <tt class="xref docutils literal"><span class="pre">'L'</span></tt>, only the lower triangular part of <tt class="docutils literal"><span class="pre">A</span></tt>
is accessed and the upper triangular part is ignored.
If <tt class="docutils literal"><span class="pre">uplo</span></tt> is <tt class="xref docutils literal"><span class="pre">'U'</span></tt>, only the upper triangular part of <tt class="docutils literal"><span class="pre">A</span></tt>
is accessed and the lower triangular part is ignored.</p>
<p>If <tt class="docutils literal"><span class="pre">p</span></tt> is not provided, a reordering from the AMD library is used.</p>
<p>The symbolic factorization is returned as an opaque C object that
can be passed to <a title="cvxopt.cholmod.numeric" class="reference internal" href="#cvxopt.cholmod.numeric"><tt class="xref docutils literal"><span class="pre">numeric</span></tt></a>.</p>
</dd></dl>
<dl class="function">
<dt id="cvxopt.cholmod.numeric">
<tt class="descclassname">cvxopt.cholmod.</tt><tt class="descname">numeric</tt><big>(</big><em>A</em>, <em>F</em><big>)</big><a class="headerlink" href="#cvxopt.cholmod.numeric" title="Permalink to this definition">¶</a></dt>
<dd><p>Performs a numeric factorization of a sparse symmetric matrix
as <a href="#equation-e-chol-ll">(3)</a> or <a href="#equation-e-chol-ldl">(4)</a>. The argument <tt class="docutils literal"><span class="pre">F</span></tt> is the
symbolic factorization computed by
<a title="cvxopt.cholmod.symbolic" class="reference internal" href="#cvxopt.cholmod.symbolic"><tt class="xref docutils literal"><span class="pre">symbolic</span></tt></a> applied to
the matrix <tt class="docutils literal"><span class="pre">A</span></tt>, or to another sparse matrix with the same sparsity
pattern and typecode, or by
<a title="cvxopt.cholmod.numeric" class="reference internal" href="#cvxopt.cholmod.numeric"><tt class="xref docutils literal"><span class="pre">numeric</span></tt></a> applied to a matrix
with the same sparsity pattern and typecode as <tt class="docutils literal"><span class="pre">A</span></tt>.</p>
<p>If <tt class="docutils literal"><span class="pre">F</span></tt> was created by a
<a title="cvxopt.cholmod.symbolic" class="reference internal" href="#cvxopt.cholmod.symbolic"><tt class="xref docutils literal"><span class="pre">symbolic</span></tt></a> with <tt class="docutils literal"><span class="pre">uplo</span></tt>
equal
to <tt class="xref docutils literal"><span class="pre">'L'</span></tt>, then only the lower triangular part of <tt class="docutils literal"><span class="pre">A</span></tt> is
accessed and the upper triangular part is ignored. If it was created
with <tt class="docutils literal"><span class="pre">uplo</span></tt> equal to <tt class="xref docutils literal"><span class="pre">'U'</span></tt>, then only the upper triangular
part of <tt class="docutils literal"><span class="pre">A</span></tt> is accessed and the lower triangular part is ignored.</p>
<p>On successful exit, the factorization is stored in <tt class="docutils literal"><span class="pre">F</span></tt>.
Raises an <tt class="xref docutils literal"><span class="pre">ArithmeticError</span></tt> if the factorization does not exist.</p>
</dd></dl>
<dl class="function">
<dt id="cvxopt.cholmod.solve">
<tt class="descclassname">cvxopt.cholmod.</tt><tt class="descname">solve</tt><big>(</big><em>F</em>, <em>B</em><span class="optional">[</span>, <em>sys = 0</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.cholmod.solve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves one of the following linear equations where <tt class="docutils literal"><span class="pre">B</span></tt> is a dense
matrix and <tt class="docutils literal"><span class="pre">F</span></tt> is the numeric factorization <a href="#equation-e-chol-ll">(3)</a>
or <a href="#equation-e-chol-ldl">(4)</a> computed by
<a title="cvxopt.cholmod.numeric" class="reference internal" href="#cvxopt.cholmod.numeric"><tt class="xref docutils literal"><span class="pre">numeric</span></tt></a>.
<tt class="docutils literal"><span class="pre">sys</span></tt> is an integer with values between 0 and 8.</p>
<table border="1" class="docutils">
<colgroup>
<col width="31%" />
<col width="69%" />
</colgroup>
<tbody valign="top">
<tr><td><tt class="docutils literal"><span class="pre">sys</span></tt></td>
<td>equation</td>
</tr>
<tr><td>0</td>
<td><img class="math" src="_images/math/6ac8905165ee8f413ccb8e724b7d8588fdfb9d1b.png" alt="AX = B"/></td>
</tr>
<tr><td>1</td>
<td><img class="math" src="_images/math/859a58c9213f8f188bab4896071aa51b0a79aeff.png" alt="LDL^TX = B"/></td>
</tr>
<tr><td>2</td>
<td><img class="math" src="_images/math/3956aa6f1545292f7bc387edc45ab2dfb06838dc.png" alt="LDX = B"/></td>
</tr>
<tr><td>3</td>
<td><img class="math" src="_images/math/937da576935afea26df6a1d4f8aeb1f2fca33ff8.png" alt="DL^TX=B"/></td>
</tr>
<tr><td>4</td>
<td><img class="math" src="_images/math/f0cd998505ca0d3bf97f0fcee4eee3e27b643072.png" alt="LX=B"/></td>
</tr>
<tr><td>5</td>
<td><img class="math" src="_images/math/0389113401db2debbec2ccd2c5ccdb3f6d80141a.png" alt="L^TX=B"/></td>
</tr>
<tr><td>6</td>
<td><img class="math" src="_images/math/0b510f9a977383df1b4d673133bceea168d95838.png" alt="DX=B"/></td>
</tr>
<tr><td>7</td>
<td><img class="math" src="_images/math/483a339f09e46bd496504b88bb3638c81fc76d1a.png" alt="P^TX=B"/></td>
</tr>
<tr><td>8</td>
<td><img class="math" src="_images/math/c48b6a571ae86596d7e87fd5b02365026c01d117.png" alt="PX=B"/></td>
</tr>
</tbody>
</table>
<p>(If <tt class="docutils literal"><span class="pre">F</span></tt> is a Cholesky factorization of the form <a href="#equation-e-chol-ll">(3)</a>,
<img class="math" src="_images/math/9ffb448918db29f2a72f8f87f421b3b3cad18f95.png" alt="D"/> is an identity matrix in this table. If <tt class="docutils literal"><span class="pre">A</span></tt> is complex,
<img class="math" src="_images/math/b623225415e8c63aaa50b34a67298cae4e2df999.png" alt="L^T"/> should be replaced by <img class="math" src="_images/math/793d926676639e8cded1205546b57f6034d84bcd.png" alt="L^H"/>.)</p>
<p>The matrix <tt class="docutils literal"><span class="pre">B</span></tt> is a dense <tt class="xref docutils literal"><span class="pre">'d'</span></tt> or <tt class="xref docutils literal"><span class="pre">'z'</span></tt> matrix, with
the same type as <tt class="docutils literal"><span class="pre">A</span></tt>. On exit it is overwritten by the solution.</p>
</dd></dl>
<dl class="function">
<dt id="cvxopt.cholmod.spsolve">
<tt class="descclassname">cvxopt.cholmod.</tt><tt class="descname">spsolve</tt><big>(</big><em>F</em>, <em>B</em><span class="optional">[</span>, <em>sys = 0</em><span class="optional">]</span><big>)</big><a class="headerlink" href="#cvxopt.cholmod.spsolve" title="Permalink to this definition">¶</a></dt>
<dd>Similar to <a title="cvxopt.cholmod.solve" class="reference internal" href="#cvxopt.cholmod.solve"><tt class="xref docutils literal"><span class="pre">solve</span></tt></a>, except that <tt class="docutils literal"><span class="pre">B</span></tt> is
a class:<cite>spmatrix</cite>, and the solution is returned as an output argument
(as an <tt class="xref docutils literal"><span class="pre">spmatrix</span></tt>). <tt class="docutils literal"><span class="pre">B</span></tt> must have the same typecode as <tt class="docutils literal"><span class="pre">A</span></tt>.</dd></dl>
<p>For the same example as above:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">X</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">8</span><span class="p">),</span> <span class="p">(</span><span class="mi">4</span><span class="p">,</span><span class="mi">2</span><span class="p">),</span> <span class="s">'d'</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">cholmod</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">F</span><span class="p">,</span> <span class="n">X</span><span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">X</span>
<span class="go">[-1.46e-01 4.88e-02]</span>
<span class="go">[ 1.33e+00 4.00e+00]</span>
<span class="go">[ 4.88e-01 1.17e+00]</span>
<span class="go">[ 2.83e+00 7.50e+00]</span>
</pre></div>
</div>
<dl class="function">
<dt id="cvxopt.cholmod.diag">
<tt class="descclassname">cvxopt.cholmod.</tt><tt class="descname">diag</tt><big>(</big><em>F</em><big>)</big><a class="headerlink" href="#cvxopt.cholmod.diag" title="Permalink to this definition">¶</a></dt>
<dd>Returns the diagonal elements of the Cholesky factor <img class="math" src="_images/math/859ccf4cd60c7bc6b8fa1afc9a42dc811a826d6f.png" alt="L"/>
in <a href="#equation-e-chol-ll">(3)</a>, as a dense matrix of the same type as <tt class="docutils literal"><span class="pre">A</span></tt>.
Note that this only applies to Cholesky factorizations. The matrix
<img class="math" src="_images/math/9ffb448918db29f2a72f8f87f421b3b3cad18f95.png" alt="D"/> in an <span class="raw-html">LDL<sup><small>T</small></sup></span>
factorization can be retrieved via <a title="cvxopt.cholmod.solve" class="reference internal" href="#cvxopt.cholmod.solve"><tt class="xref docutils literal"><span class="pre">solve</span></tt></a>
with <tt class="docutils literal"><span class="pre">sys</span></tt> equal to 6.</dd></dl>
<p>In the functions listed above, the default values of the control
parameters described in the CHOLMOD user guide are used, except for
<tt class="xref docutils literal"><span class="pre">Common->print</span></tt> which is set to 0 instead of 3 and
<tt class="xref docutils literal"><span class="pre">Common->supernodal</span></tt> which is set to 2 instead of 1.
These parameters (and a few others) can be modified by making an
entry in the dictionary <tt class="xref docutils literal"><span class="pre">cholmod.options</span></tt>.
The meaning of these parameters is as follows.</p>
<dl class="docutils">
<dt><tt class="xref docutils literal"><span class="pre">options['supernodal']</span></tt></dt>
<dd>If equal to 0, a factorization <a href="#equation-e-chol-ldl">(4)</a> is computed using a
simplicial algorithm. If equal to 2, a factorization <a href="#equation-e-chol-ll">(3)</a>
is computed using a supernodal algorithm. If equal to 1, the most
efficient of the two factorizations is selected, based on the sparsity
pattern. Default: 2.</dd>
<dt><tt class="xref docutils literal"><span class="pre">options['print']</span></tt></dt>
<dd>A nonnegative integer that controls the amount of output printed to
the screen. Default: 0 (no output).</dd>
</dl>
<p>As an example that illustrates <a title="cvxopt.cholmod.diag" class="reference internal" href="#cvxopt.cholmod.diag"><tt class="xref docutils literal"><span class="pre">diag</span></tt></a> and the
use of <tt class="xref docutils literal"><span class="pre">cholmod.options</span></tt>, we compute the logarithm of the determinant
of the coefficient matrix in <a href="#equation-e-A-pd">(2)</a> by two methods.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">import</span> <span class="nn">math</span>
<span class="gp">>>> </span><span class="kn">from</span> <span class="nn">cvxopt.cholmod</span> <span class="kn">import</span> <span class="n">options</span>
<span class="gp">>>> </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">log</span>
<span class="gp">>>> </span><span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="mf">2.0</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">cholmod</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">F</span><span class="p">)))</span>
<span class="go">5.50533153593</span>
<span class="gp">>>> </span><span class="n">options</span><span class="p">[</span><span class="s">'supernodal'</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span>
<span class="gp">>>> </span><span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">Di</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">4</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
<span class="gp">>>> </span><span class="n">cholmod</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">F</span><span class="p">,</span> <span class="n">Di</span><span class="p">,</span> <span class="n">sys</span><span class="o">=</span><span class="mi">6</span><span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="o">-</span><span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">Di</span><span class="p">))</span>
<span class="go">5.50533153593</span>
</pre></div>
</div>
</div>
<div class="section" id="example-covariance-selection">
<h2>Example: Covariance Selection<a class="headerlink" href="#example-covariance-selection" title="Permalink to this headline">¶</a></h2>
<p>This example illustrates the use of the routines for sparse Cholesky
factorization. We consider the problem</p>
<div class="math" id="equation-e-covsel">
<p><span class="eqno">(5)</span><img src="_images/math/9a52c242e18c58867c799409c61df65280145c29.png" alt="\newcommand{\Tr}{\mathop{\bf tr}}
\begin{array}{ll}
\mbox{minimize} & -\log\det K + \Tr(KY) \\
\mbox{subject to} & K_{ij}=0,\quad (i,j) \not \in S.
\end{array}" /></p>
</div><p>The optimization variable is a symmetric matrix <img class="math" src="_images/math/dfb064112b6c94470339f6571f69d07afc1c024c.png" alt="K"/> of order
<img class="math" src="_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/> and the domain of the problem is the set of positive definite
matrices. The matrix <img class="math" src="_images/math/ce58e4af225c93d08606c26554caaa5ae32edeba.png" alt="Y"/> and the index set <img class="math" src="_images/math/ad28c83c99a8fd0dd2e2e594c9d02ee532765a0a.png" alt="S"/> are given.
We assume that all the diagonal positions are included in <img class="math" src="_images/math/ad28c83c99a8fd0dd2e2e594c9d02ee532765a0a.png" alt="S"/>.
This problem arises in maximum likelihood estimation of the covariance
matrix of a zero-mean normal distribution, with constraints
that specify that pairs of variables are conditionally independent.</p>
<p>We can express <img class="math" src="_images/math/dfb064112b6c94470339f6571f69d07afc1c024c.png" alt="K"/> as</p>
<div class="math">
<p><img src="_images/math/67b73ba16f26452982e53f9f5da5d8c246cbd4de.png" alt="\newcommand{\diag}{\mathop{\bf diag}}
K(x) = E_1\diag(x)E_2^T+E_2\diag(x)E_1^T" /></p>
</div><p>where <img class="math" src="_images/math/26eeb5258ca5099acf8fe96b2a1049c48c89a5e6.png" alt="x"/> are the nonzero elements in the lower triangular part of
<img class="math" src="_images/math/dfb064112b6c94470339f6571f69d07afc1c024c.png" alt="K"/>, with the diagonal elements scaled by 1/2, and</p>
<div class="math">
<p><img src="_images/math/77fbceca18518dc4b8b4994c8c5ce16d02d68a4c.png" alt="E_1 = \left[ \begin{array}{cccc}
e_{i_1} & e_{i_2} & \cdots & e_{i_q} \end{array}\right], \qquad
E_2 = \left[ \begin{array}{cccc}
e_{j_1} & e_{j_2} & \cdots & e_{j_q} \end{array}\right]," /></p>
</div><p>where (<img class="math" src="_images/math/d918c0ad087769bca058188217bdc29e865d685c.png" alt="i_k"/>, <img class="math" src="_images/math/5be7f9688416e48fd04f912dd6e1dff4749f9db0.png" alt="j_k"/>) are the positions of the nonzero entries
in the lower-triangular part of <img class="math" src="_images/math/dfb064112b6c94470339f6571f69d07afc1c024c.png" alt="K"/>. With this notation, we can
solve problem <a href="#equation-e-covsel">(5)</a> by solving the unconstrained problem</p>
<div class="math">
<p><img src="_images/math/73ef4db286e93bab211267c97ffde30973795f5a.png" alt="\newcommand{\Tr}{\mathop{\bf tr}}
\begin{array}{ll}
\mbox{minimize} & f(x) = -\log\det K(x) + \Tr(K(x)Y).
\end{array}" /></p>
</div><p>The code below implements Newton’s method with a backtracking line search.
The gradient and Hessian of the objective function are given by</p>
<div class="math">
<p><img src="_images/math/3437890f06999b22bb720f0cbb50ad44cc671af1.png" alt="\newcommand{\diag}{\mathop{\bf diag}}
\begin{split}
\nabla f(x)
& = 2 \diag( E_1^T (Y - K(x)^{-1}) E_2)) \\
& = 2\diag(Y_{IJ} - \left(K(x)^{-1}\right)_{IJ}) \\
\nabla^2 f(x)
& = 2 (E_1^T K(x)^{-1} E_1) \circ (E_2^T K(x)^{-1} E_2)
+ 2 (E_1^T K(x)^{-1} E_2) \circ (E_2^T K(x)^{-1} E_1) \\
& = 2 \left(K(x)^{-1}\right)_{II} \circ \left(K(x)^{-1}\right)_{JJ}
+ 2 \left(K(x)^{-1}\right)_{IJ} \circ
\left(K(x)^{-1}\right)_{JI},
\end{split}" /></p>
</div><p>where <img class="math" src="_images/math/7f3b4ea3b8abe9ced6095c7d75f207a1127ee843.png" alt="\circ"/> denotes Hadamard product.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">log</span><span class="p">,</span> <span class="n">mul</span><span class="p">,</span> <span class="n">blas</span><span class="p">,</span> <span class="n">lapack</span><span class="p">,</span> <span class="n">amd</span><span class="p">,</span> <span class="n">cholmod</span>
<span class="k">def</span> <span class="nf">covsel</span><span class="p">(</span><span class="n">Y</span><span class="p">):</span>
<span class="sd">"""</span>
<span class="sd"> Returns the solution of</span>
<span class="sd"> minimize -logdet K + Tr(KY)</span>
<span class="sd"> subject to K_{ij}=0, (i,j) not in indices listed in I,J.</span>
<span class="sd"> Y is a symmetric sparse matrix with nonzero diagonal elements.</span>
<span class="sd"> I = Y.I, J = Y.J.</span>
<span class="sd"> """</span>
<span class="n">I</span><span class="p">,</span> <span class="n">J</span> <span class="o">=</span> <span class="n">Y</span><span class="o">.</span><span class="n">I</span><span class="p">,</span> <span class="n">Y</span><span class="o">.</span><span class="n">J</span>
<span class="n">n</span><span class="p">,</span> <span class="n">m</span> <span class="o">=</span> <span class="n">Y</span><span class="o">.</span><span class="n">size</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="nb">len</span><span class="p">(</span><span class="n">I</span><span class="p">)</span>
<span class="n">N</span> <span class="o">=</span> <span class="n">I</span> <span class="o">+</span> <span class="n">J</span><span class="o">*</span><span class="n">n</span> <span class="c"># non-zero positions for one-argument indexing</span>
<span class="n">D</span> <span class="o">=</span> <span class="p">[</span><span class="n">k</span> <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">m</span><span class="p">)</span> <span class="k">if</span> <span class="n">I</span><span class="p">[</span><span class="n">k</span><span class="p">]</span><span class="o">==</span><span class="n">J</span><span class="p">[</span><span class="n">k</span><span class="p">]]</span> <span class="c"># position of diagonal elements</span>
<span class="c"># starting point: symmetric identity with nonzero pattern I,J</span>
<span class="n">K</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>
<span class="n">K</span><span class="p">[::</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mi">1</span>
<span class="c"># Kn is used in the line search</span>
<span class="n">Kn</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>
<span class="c"># symbolic factorization of K</span>
<span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">K</span><span class="p">)</span>
<span class="c"># Kinv will be the inverse of K</span>
<span class="n">Kinv</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">))</span>
<span class="k">for</span> <span class="n">iters</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="mi">100</span><span class="p">):</span>
<span class="c"># numeric factorization of K</span>
<span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">K</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="n">d</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">F</span><span class="p">)</span>
<span class="c"># compute Kinv by solving K*X = I</span>
<span class="n">Kinv</span><span class="p">[:]</span> <span class="o">=</span> <span class="mi">0</span>
<span class="n">Kinv</span><span class="p">[::</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mi">1</span>
<span class="n">cholmod</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">F</span><span class="p">,</span> <span class="n">Kinv</span><span class="p">)</span>
<span class="c"># solve Newton system</span>
<span class="n">grad</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="p">(</span><span class="n">Y</span><span class="o">.</span><span class="n">V</span> <span class="o">-</span> <span class="n">Kinv</span><span class="p">[</span><span class="n">N</span><span class="p">])</span>
<span class="n">hess</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="p">(</span><span class="n">mul</span><span class="p">(</span><span class="n">Kinv</span><span class="p">[</span><span class="n">I</span><span class="p">,</span><span class="n">J</span><span class="p">],</span><span class="n">Kinv</span><span class="p">[</span><span class="n">J</span><span class="p">,</span><span class="n">I</span><span class="p">])</span> <span class="o">+</span> <span class="n">mul</span><span class="p">(</span><span class="n">Kinv</span><span class="p">[</span><span class="n">I</span><span class="p">,</span><span class="n">I</span><span class="p">],</span><span class="n">Kinv</span><span class="p">[</span><span class="n">J</span><span class="p">,</span><span class="n">J</span><span class="p">]))</span>
<span class="n">v</span> <span class="o">=</span> <span class="o">-</span><span class="n">grad</span>
<span class="n">lapack</span><span class="o">.</span><span class="n">posv</span><span class="p">(</span><span class="n">hess</span><span class="p">,</span><span class="n">v</span><span class="p">)</span>
<span class="c"># stopping criterion</span>
<span class="n">sqntdecr</span> <span class="o">=</span> <span class="o">-</span><span class="n">blas</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">grad</span><span class="p">,</span><span class="n">v</span><span class="p">)</span>
<span class="k">print</span> <span class="s">"Newton decrement squared:</span><span class="si">%- 7.5e</span><span class="s">"</span> <span class="o">%</span><span class="n">sqntdecr</span>
<span class="k">if</span> <span class="p">(</span><span class="n">sqntdecr</span> <span class="o"><</span> <span class="mf">1e-12</span><span class="p">):</span>
<span class="k">print</span> <span class="s">"number of iterations: "</span><span class="p">,</span> <span class="n">iters</span><span class="o">+</span><span class="mi">1</span>
<span class="k">break</span>
<span class="c"># line search</span>
<span class="n">dx</span> <span class="o">=</span> <span class="o">+</span><span class="n">v</span>
<span class="n">dx</span><span class="p">[</span><span class="n">D</span><span class="p">]</span> <span class="o">*=</span> <span class="mi">2</span> <span class="c"># scale the diagonal elems</span>
<span class="n">f</span> <span class="o">=</span> <span class="o">-</span><span class="mf">2.0</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">))</span> <span class="c"># f = -log det K</span>
<span class="n">s</span> <span class="o">=</span> <span class="mi">1</span>
<span class="k">for</span> <span class="n">lsiter</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="mi">50</span><span class="p">):</span>
<span class="n">Kn</span><span class="o">.</span><span class="n">V</span> <span class="o">=</span> <span class="n">K</span><span class="o">.</span><span class="n">V</span> <span class="o">+</span> <span class="n">s</span><span class="o">*</span><span class="n">dx</span>
<span class="k">try</span><span class="p">:</span>
<span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">Kn</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="k">except</span> <span class="ne">ArithmeticError</span><span class="p">:</span>
<span class="n">s</span> <span class="o">*=</span> <span class="mf">0.5</span>
<span class="k">else</span><span class="p">:</span>
<span class="n">d</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">F</span><span class="p">)</span>
<span class="n">fn</span> <span class="o">=</span> <span class="o">-</span><span class="mf">2.0</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">))</span> <span class="o">+</span> <span class="mi">2</span><span class="o">*</span><span class="n">s</span><span class="o">*</span><span class="n">blas</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">v</span><span class="p">,</span><span class="n">Y</span><span class="o">.</span><span class="n">V</span><span class="p">)</span>
<span class="k">if</span> <span class="p">(</span><span class="n">fn</span> <span class="o"><</span> <span class="n">f</span> <span class="o">-</span> <span class="mf">0.01</span><span class="o">*</span><span class="n">s</span><span class="o">*</span><span class="n">sqntdecr</span><span class="p">):</span>
<span class="k">break</span>
<span class="n">s</span> <span class="o">*=</span> <span class="mf">0.5</span>
<span class="n">K</span><span class="o">.</span><span class="n">V</span> <span class="o">=</span> <span class="n">Kn</span><span class="o">.</span><span class="n">V</span>
<span class="k">return</span> <span class="n">K</span>
</pre></div>
</div>
</div>
</div>
</div>
</div>
</div>
<div class="sphinxsidebar">
<div class="sphinxsidebarwrapper">
<h3><a href="index.html">Table Of Contents</a></h3>
<ul>
<li><a class="reference external" href="#">Sparse Linear Equations</a><ul>
<li><a class="reference external" href="#matrix-orderings">Matrix Orderings</a></li>
<li><a class="reference external" href="#general-linear-equations">General Linear Equations</a></li>
<li><a class="reference external" href="#positive-definite-linear-equations">Positive Definite Linear Equations</a></li>
<li><a class="reference external" href="#example-covariance-selection">Example: Covariance Selection</a></li>
</ul>
</li>
</ul>
<h4>Previous topic</h4>
<p class="topless"><a href="fftw.html"
title="previous chapter">Discrete Transforms</a></p>
<h4>Next topic</h4>
<p class="topless"><a href="coneprog.html"
title="next chapter">Cone Programming</a></p>
<div id="searchbox" style="display: none">
<h3>Quick search</h3>
<form class="search" action="search.html" method="get">
<input type="text" name="q" size="18" />
<input type="submit" value="Go" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
<p class="searchtip" style="font-size: 90%">
Enter search terms or a module, class or function name.
</p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="coneprog.html" title="Cone Programming"
>next</a></li>
<li class="right" >
<a href="fftw.html" title="Discrete Transforms"
>previous</a> |</li>
<li><a href="http://abel.ee.ucla.edu/cvxopt">CVXOPT home</a> |</li>
<li><a href="index.html">user's guide</a> </li>
</ul>
</div>
<div class="footer">
Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 0.6.5.
</div>
</body>
</html>
