Changeset 901

Show
Ignore:
Timestamp:
09/23/05 11:13:49 (3 years ago)
Author:
rkern
Message:

r3086@Blasphemy: kern | 2005-09-21 21:49:31 -0700
more scipy_tutorial fixes

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • nbdoc/trunk/scipy_tutorial/chapter01.pybk

    r845 r901  
    2626the command 
    2727</para> 
    28 <para>XXX: ipython 
     28<programlisting>XXX: ipython 
    2929>>> from scipy import * 
    30 </para> 
     30</programlisting> 
    3131 
    3232<section> 
    …… 
    5151in that module is printed. For example: 
    5252</para> 
    53 <para>XXX: ipython 
     53<programlisting>XXX: ipython 
    5454>>> info(optimize.fmin) 
    5555 fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, 
    …… 
    8888  printmessg -- non-zero to print convergence messages. 
    8989   
    90 </para> 
     90</programlisting> 
    9191<para>Another useful command is <code>source</code>. When given a function written in 
    9292Python as an argument, it prints out a listing of the source code for that 

  • nbdoc/trunk/scipy_tutorial/chapter02.pybk

    r854 r901  
    3131have been modified to return complex numbers instead of NaN's where appropriate. 
    3232</para> 
    33 <para>XXX: ipython 
    34 >>> scipy.sqrt(-1) </para> 
     33<programlisting>XXX: ipython 
     34>>> scipy.sqrt(-1) </programlisting> 
    3535</section> 
    3636<section> 
    …… 
    122122its usage:  
    123123</para> 
    124 <para>XXX:  
     124<programlisting>XXX:  
    125125&gt;&gt;&gt; mgrid[0:5,0:5] 
    126126array([[[0, 0, 0, 0, 0], 
    …… 
    143143        [ 0.    ,  1.6667,  3.3333,  5.    ], 
    144144        [ 0.    ,  1.6667,  3.3333,  5.    ]]]) 
    145 </para> 
     145</programlisting> 
    146146<para>Having meshed arrays like this is sometimes very useful. However, it is 
    147147not always needed just to evaluate some N-dimensional function over a grid due 
    …… 
    179179differentiated, and evaluated. It even prints like a polynomial: 
    180180</para> 
    181 <para> XXX: ipython 
     181<programlisting> XXX: ipython 
    182182&gt;&gt;&gt; p = 
    183183poly1d([3,4,5]) 
    …… 
    195195&gt;&gt;&gt; p([4,5]) 
    196196array([ 69, 100]) 
    197 </para> 
     197</programlisting> 
    198198<para>The other way to handle polynomials is as an array of coefficients with 
    199199the first element of the array giving the coefficient of the highest power. 
    …… 
    212212have a Python function named addsubtract defined as: 
    213213</para> 
    214 <para>XXX: ipython 
     214<programlisting>XXX: ipython 
    215215&gt;&gt;&gt; def addsubtract(a,b): 
    216216    if a &gt; b: 
    …… 
    218218    else: 
    219219        return a + b 
    220 </para> 
     220</programlisting> 
    221221<para>which defines a function of two scalar variables and returns a scalar 
    222222result. The class vectorize can be used to "vectorize" this function so that 
    223223returns a function which takes array arguments and returns an array result: 
    224224</para> 
    225 <para> XXX: ipython 
     225<programlisting> XXX: ipython 
    226226&gt;&gt;&gt; vec_addsubtract([0,3,6,9],[1,3,5,7]) 
    227227array([1, 6, 1, 2]) 
    228 </para> 
     228</programlisting> 
    229229<para>This particular function could have been written in vector form without 
    230230the use of vectorize. But, what if the function you have written is the result 
    …… 
    254254corresponding to the first condition in condlist that is true. For example  
    255255</para> 
    256 <para> XXX: ipython 
     256<programlisting> XXX: ipython 
    257257&gt;&gt;&gt; x = r_[-2:3] 
    258258&gt;&gt;&gt; x 
    …… 
    260260&gt;&gt;&gt; select([x &gt; 3, x &gt;= 0],[0,x+2]) 
    261261array([0, 0, 2, 3, 4]) 
    262 </para> 
     262</programlisting> 
    263263</section> 
    264264</section> 

  • nbdoc/trunk/scipy_tutorial/chapter04.pybk

    r853 r901  
    55module is provided by the help command: 
    66</para> 
    7 <para>XXX: ipython 
     7<programlisting>XXX: ipython 
    88&gt;&gt;&gt; help(integrate) 
    99Methods for Integrating Functions 
    …… 
    1818  See the orthogonal module (integrate.orthogonal) for Gaussian 
    1919     quadrature roots and weights. 
    20 </para> 
     20</programlisting> 
    2121<section> 
    2222<title>General integration (integrate.quad)</title> 
    …… 
    3030This could be computed using quad: 
    3131</para> 
    32 <para>XXX: ipython 
     32<programlisting>XXX: ipython 
    3333&gt;&gt;&gt; result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5) 
    3434&gt;&gt;&gt; print result 
    …… 
    4242&gt;&gt;&gt; print abs(result[0]-I) 
    43431.03761443881e-11  
    44 </para> 
     44</programlisting> 
    4545<para>The first argument to quad is a "callable" Python object (i.e a function, 
    4646method, or class instance). Notice the use of a lambda-function in this case as 
    …… 
    7272defining a new function vec_expint based on the routine quad:  
    7373</para> 
    74 <para> XXX: ipython 
     74<programlisting> XXX: ipython 
    7575&gt;&gt;&gt; from integrate import quad, Inf 
    7676&gt;&gt;&gt; def integrand(t,n,x): 
    …… 
    8686&gt;&gt;&gt; special.expn(3,arange(1.0,4.0,0.5)) 
    8787array([ 0.1097,  0.0567,  0.0301,  0.0163,  0.0089,  0.0049]) 
    88 </para> 
     88</programlisting> 
    8989<para>The function which is integrated can even use the quad argument (though 
    9090the error bound may underestimate the error due to possible numerical error in 
    …… 
    9494dx=\frac{1}{n}</ipython-equation>. 
    9595</para> 
    96 <para>XXX: ipython 
     96<programlisting>XXX: ipython 
    9797&gt;&gt;&gt; result = quad(lambda x: expint(3, x), 0, Inf) 
    9898&gt;&gt;&gt; print result 
    …… 
    105105&gt;&gt;&gt; print I3 - result[0] 
    1061068.77306560731e-11      
    107 </para> 
     107</programlisting> 
    108108<para>This last example shows that multiple integration can be handled using 
    109109repeated calls to quad. The mechanics of this for double and triple integration 
    …… 
    115115<ipython-inlineequation>I_{n}</ipython-inlineequation> is shown below: 
    116116</para> 
    117 <para>XXX: ipython 
     117<programlisting>XXX: ipython 
    118118&gt;&gt;&gt; from __future__ import nested_scopes 
    119119&gt;&gt;&gt; from integrate import quad, dblquad, Inf 
    …… 
    127127&gt;&gt;&gt; print I(2) 
    128128(0.49999999999857514, 1.8855523253868967e-09) 
    129 </para> 
     129</programlisting> 
    130130</section> 
    131131<section> 
    …… 
    237237The following example illustrates the use of odeint including the usage of the 
    238238Dfun option which allows the user to specify a gradient (with respect to 
    239 <ipython-inlineequation>\mathbf{y})</ipython-inlineequation> of the function, 
     239<ipython-inlineequation>\mathbf{y}</ipython-inlineequation>) of the function, 
    240240<ipython-inlineequation>\mathbf{f}\left(\mathbf{y},t\right)</ipython-inlineequation>. 
    241241</para> 
    242 <para>XXX: ipython 
     242<programlisting>XXX: ipython 
    243243&gt;&gt;&gt; from integrate import odeint 
    244244&gt;&gt;&gt; from special import gamma, airy 
    …… 
    268268&gt;&gt;&gt; print y2[:36:6,1] 
    269269[ 0.355028  0.339511  0.324067  0.308763  0.293658  0.278806] 
    270 </para> 
     270</programlisting> 
    271271 
    272272</section> 

  • nbdoc/trunk/scipy_tutorial/chapter05.pybk

    r845 r901  
    55pydoc.help): 
    66</para> 
    7 <para>XXX: ipython 
     7<programlisting>XXX: ipython 
    88&gt;&gt;&gt; info(optimize) 
    99 Optimization Tools 
    …… 
    4242 
    4343  fixed_point -- Single-variable fixed-point solver. 
    44 </para> 
     44</programlisting> 
    4545<para> 
    4646 The first four algorithms are unconstrained minimization algorithms (fmin: 
    …… 
    6969minimum can be found using the fmin routine as shown in the example below: 
    7070</para> 
    71 <para>XXX: ipython 
     71<programlisting>XXX: ipython 
    7272&gt;&gt;&gt; from scipy.optimize import fmin 
    7373&gt;&gt;&gt; def rosen(x):  # The Rosenbrock function 
    …… 
    8383&gt;&gt;&gt; print xopt 
    8484[ 1.  1.  1.  1.  1.] 
    85 </para> 
     85</programlisting> 
    8686<para>Another optimization algorithm that needs only function calls to find the 
    8787minimum is Powell's method available as optimize.fmin_powell. 
    …… 
    123123gradient is constructed by the code-segment: 
    124124</para> 
    125 <para>XXX: ipython 
     125<programlisting>XXX: ipython 
    126126&gt;&gt;&gt; def rosen_der(x): 
    127127        xm = x[1:-1] 
    …… 
    133133        der[-1] = 200*(x[-1]-x[-2]**2) 
    134134        return der 
    135 </para> 
     135</programlisting> 
    136136<para> 
    137137The calling signature for the BFGS minimization algorithm is similar to fmin 
    …… 
    139139in the following example which minimizes the Rosenbrock function. 
    140140</para> 
    141 <para>XXX: ipython 
     141<programlisting>XXX: ipython 
    142142&gt;&gt;&gt; from scipy.optimize import fmin_bfgs 
    143143 
    …… 
    151151&gt;&gt;&gt; print xopt 
    152152[ 1.  1.  1.  1.  1.]  
    153 </para> 
     153</programlisting> 
    154154</section> 
    155155 
    …… 
    217217the following example: 
    218218</para> 
    219 <para>XXX: ipython 
     219<programlisting>XXX: ipython 
    220220&gt;&gt;&gt; from scipy.optimize import fmin_ncg 
    221221&gt;&gt;&gt; def rosen_hess(x): 
    …… 
    239239&gt;&gt;&gt; print xopt 
    240240[ 0.9999  0.9999  0.9998  0.9996  0.9991] 
    241 </para> 
     241</programlisting> 
    242242</section> 
    243243<section> 
    …… 
    271271fhess_p keyword to minimize the Rosenbrock function using fmin_ncg follows: 
    272272</para> 
    273 <para>XXX: ipython 
     273<programlisting>XXX: ipython 
    274274&gt;&gt;&gt; from scipy.optimize import fmin_ncg 
    275275&gt;&gt;&gt; def rosen_hess_p(x,p): 
    …… 
    292292&gt;&gt;&gt; print xopt 
    293293[ 1.      1.      1.      0.9999  0.9999] 
    294 </para> 
     294</programlisting> 
    295295</section> 
    296296</section> 
    …… 
    350350example and a plot of the results is shown in Figure [fig:least_squares_fit]. 
    351351</para> 
    352 <para>XXX: ipython and figure 
     352<programlisting>XXX: ipython and figure 
    353353&gt;&gt;&gt; x = arange(0,6e-2,6e-2/30) 
    354354&gt;&gt;&gt; A,k,theta = 10, 1.0/3e-2, pi/6 
    …… 
    381381&gt;&gt;&gt; legend(['Fit', 'Noisy', 'True']) 
    382382&gt;&gt;&gt; gist.eps('leastsqfit')   # Make epsi file. 
    383 </para> 
     383</programlisting> 
    384384</section> 
    385385 
    …… 
    426426<ipython-inlineequation>x_{\textrm{min}}=5.3314</ipython-inlineequation>: 
    427427</para> 
    428 <para> XXX: ipython 
     428<programlisting> XXX: ipython 
    429429&gt;&gt;&gt; from scipy.special import j1 
    430430 
    …… 
    433433&gt;&gt;&gt; print xmin 
    4344345.33144184241 
    435 </para> 
     435</programlisting> 
    436436</section> 
    437437</section> 
    …… 
    456456<ipython-inlineequation>x_{0}=6.5041,\, x_{1}=0.9084</ipython-inlineequation>. 
    457457</para> 
    458 <para>XXX: ipython 
     458<programlisting>XXX: ipython 
    459459&gt;&gt;&gt; def func(x): 
    460460        return x + 2*cos(x) 
    …… 
    473473&gt;&gt;&gt; print x02 
    474474[ 6.5041  0.9084] 
    475 </para> 
     475</programlisting> 
    476476</section> 
    477477<section> 

  • nbdoc/trunk/scipy_tutorial/chapter06.pybk

    r845 r901  
    1919be specified at instantiation time. The following example demonstrates it's use.  
    2020</para> 
    21 <para> XXX: ipython 
     21<programlisting> XXX: ipython 
    2222&gt;&gt;&gt; x = arange(0,10) 
    2323&gt;&gt;&gt; y = exp(-x/3.0) 
    …… 
    3838&gt;&gt;&gt; xnew = arange(0,9,0.1) 
    3939&gt;&gt;&gt; xplt.plot(x,y,'x',xnew,f(xnew),'.') 
    40 </para> 
     40</programlisting> 
    4141<para> 
    4242Figure shows the result: [Graphics file: inter_1d.epsi] 
    …… 
    9696[fig:spline-1d]Examples of using cubic-spline interpolation. 
    9797</para> 
    98 <para> XXX: ipython 
     98<programlisting> XXX: ipython 
    9999&gt;&gt;&gt; # Cubic-spline 
    100100&gt;&gt;&gt; x = arange(0,2*pi+pi/4,2*pi/8) 
    …… 
    149149&gt;&gt;&gt; xplt.title('Spline of parametrically-defined curve') 
    150150&gt;&gt;&gt; xplt.eps('interp_cubic_param') 
    151 </para> 
     151</programlisting> 
    152152</section> 
    153153 
    …… 
    190190indexing objects passed in mgrid[].  
    191191</para> 
    192 <para> XXX: ipython 
     192<programlisting> XXX: ipython 
    193193&gt;&gt;&gt; # Define function over sparse 20x20 grid 
    194194&gt;&gt;&gt; x,y = mgrid[-1:1:20j,-1:1:20j] 
    …… 
    205205&gt;&gt;&gt; xplt.title3("Interpolated function.") 
    206206&gt;&gt;&gt; xplt.eps("2d_interp") 
    207 </para> 
     207</programlisting> 
    208208<para> 
    209209[Graphics file: 2d_func.epsi] 

  • nbdoc/trunk/scipy_tutorial/chapter07.pybk

    r845 r901  
    9595and allows for choosing mirror-symmetric boundary conditions.  
    9696</para> 
    97 <para>XXX: ipython 
     97<programlisting>XXX: ipython 
    9898&gt;&gt;&gt; image = lena().astype(Float32) 
    9999&gt;&gt;&gt; derfilt = array([1.0,-2,1.0],Float32) 
    …… 
    112112&gt;&gt;&gt; xplt.title('Output of spline edge filter') 
    113113&gt;&gt;&gt; xplt.eps('lena_edge') 
    114 </para> 
     114</programlisting> 
    115115<para> 
    116116[Graphics file: lena_image.epsi] 

  • nbdoc/trunk/scipy_tutorial/chapter10.pybk

    r853 r901  
    5050 The following example demonstrates this computation in SciPy 
    5151</para> 
    52 <para> XXX: ipython 
     52<programlisting> XXX: ipython 
    5353&gt;&gt;&gt; A = mat('[1 3 5; 2 5 1; 2 3 8]') 
    5454&gt;&gt;&gt; A 
    …… 
    6464       [ 0.56,  0.08, -0.36], 
    6565       [ 0.16, -0.12,  0.04]]) 
    66 </para> 
     66</programlisting> 
    6767</section> 
    6868 
    …… 
    103103answer as shown in the following example:  
    104104</para> 
    105 <para>XXX: ipython 
     105<programlisting>XXX: ipython 
    106106&gt;&gt;&gt; A = mat('[1 3 5; 2 5 1; 2 3 8]') 
    107107&gt;&gt;&gt; b = mat('[10;8;3]') 
    …… 
    114114       [ 5.16], 
    115115       [ 0.76]]) 
    116 </para> 
     116</programlisting> 
    117117</section> 
    118118 
    119119<section> 
    120120<title>Finding Determinant</title> 
    121 <para>The determinant of a square matrix \mathbf{A} is often denoted 
    122 \left|\mathbf{A}\right| and is a quantity often used in linear algebra. Suppose 
    123 a_{ij} are the elements of the matrix \mathbf{A} and let 
    124 M_{ij}=\left|\mathbf{A}_{ij}\right| be the determinant of the matrix left by 
    125 removing the i^{\textrm{th}} row and j^{\textrm{th}}column from \mathbf{A}. Then 
    126 for any row i,\left|\mathbf{A}\right|=\sum_{j}\left(-1\right)^{i+j}a_{ij}M_{ij}. 
     121<para>The determinant of a square matrix 
     122<ipython-inlineequation>\mathbf{A}</ipython-inlineequation> is often denoted 
     123<ipython-inlineequation>\left|\mathbf{A}\right|</ipython-inlineequation> and is a quantity often used in linear algebra. Suppose 
     124<ipython-inlineequation>a_{ij}</ipython-inlineequation> are the elements of the matrix <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> and let 
     125<ipython-inlineequation>M_{ij}=\left|\mathbf{A}_{ij}\right|</ipython-inlineequation> be the determinant of the matrix left by 
     126removing the <ipython-inlineequation>i^{\textrm{th}}</ipython-inlineequation> row and <ipython-inlineequation>j^{\textrm{th}}column</ipython-inlineequation> from <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>. Then 
     127for any row <ipython-equation>i,\left|\mathbf{A}\right|=\sum_{j}\left(-1\right)^{i+j}a_{ij}M_{ij}</ipython-equation>. 
    127128This is a recursive way to define the determinant where the base case is defined 
    128129by accepting that the determinant of a 1\times1 matrix is the only matrix 
    129130element. In SciPy the determinant can be calculated with linalg.det. For 
    130 example, the determinant of \mathbf{A=}\left[\begin{array}{ccc} 
     131example, the determinant of  
     132<ipython-equation>\mathbf{A=}\left[\begin{array}{ccc} 
    1311331 &amp; 3 &amp; 5\\ 
    1321342 &amp; 5 &amp; 1\\ 
    133 2 &amp; 3 &amp; 8\end{array}\right]is \left|\mathbf{A}\right| In SciPy this is 
     1352 &amp; 3 &amp; 8\end{array}\right]</ipython-equation> is  
     136<ipython-equation verb="1">\begin{eqnarray*} 
     137\left|\mathbf{A}\right| &amp; = &amp; 1\left|\begin{array}{cc} 
     1385 &amp; 1\\ 
     1393 &amp; 8\end{array}\right|-3\left|\begin{array}{cc} 
     1402 &amp; 1\\ 
     1412 &amp; 8\end{array}\right|+5\left|\begin{array}{cc} 
     1422 &amp; 5\\ 
     1432 &amp; 3\end{array}\right|\\ 
     144 &amp; = &amp; 1\left(5\cdot8-3\cdot1\right)-3\left(2\cdot8-2\cdot1\right)+5\left(2\cdot3-2\cdot5\right)=-25 
     145.\end{eqnarray*} 
     146</ipython-equation> In SciPy this is 
    134147computed as shown in this example: 
    135148</para> 
    136 <para> XXX: ipython 
     149<programlisting> XXX: ipython 
    137150&gt;&gt;&gt; A = mat('[1 3 5; 2 5 1; 2 3 8]') 
    138151&gt;&gt;&gt; linalg.det(A) 
    139152-25.000000000000004 
    140 </para> 
     153</programlisting> 
    141154</section> 
    142155 
    …… 
    150163</para> 
    151164<para> 
    152 For vector x, the order parameter can be any real number including inf or -inf. The computed norm is \left\Vert \mathbf{x}\right\Vert =\left\{ \begin{array}{cc} 
     165For vector x, the order parameter can be any real number including inf or -inf. The computed norm is  
     166<ipython-equation>\left\Vert \mathbf{x}\right\Vert =\left\{ \begin{array}{cc} 
    153167\max\left|x_{i}\right| &amp; \textrm{ord}=\textrm{inf}\\ 
    154168\min\left|x_{i}\right| &amp; \textrm{ord}=-\textrm{inf}\\ 
    155169\left(\sum_{i}\left|x_{i}\right|^{\textrm{ord}}\right)^{1/\textrm{ord}} &amp; 
    156 \left|\textrm{ord}\right|&lt;\infty.\end{array}\right. 
     170\left|\textrm{ord}\right|&lt;\infty.\end{array}\right.</ipython-equation>  
    157171</para> 
    158172<para> 
    159 For matrix \mathbf{A} the only valid values for norm are \pm2,\pm1, \pm inf, and 
    160 'fro' (or 'f') Thus, \left\Vert \mathbf{A}\right\Vert =\left\{ \begin{array}{cc} 
     173For matrix <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> the only valid values for norm are <ipython-inlineequation>\pm2</ipython-inlineequation>,<ipython-inlineequation>\pm1</ipython-inlineequation>, <ipython-inlineequation>\pm inf</ipython-inlineequation>, and 
     174'fro' (or 'f'). Thus, <ipython-equation>\left\Vert \mathbf{A}\right\Vert =\left\{ \begin{array}{cc} 
    161175\max_{i}\sum_{j}\left|a_{ij}\right| &amp; \textrm{ord}=\textrm{inf}\\ 
    162176\min_{i}\sum_{j}\left|a_{ij}\right| &amp; \textrm{ord}=-\textrm{inf}\\ 
    …… 
    166180\min\sigma_{i} &amp; \textrm{ord}=-2\\ 
    167181\sqrt{\textrm{trace}\left(\mathbf{A}^{H}\mathbf{A}\right)} &amp; 
    168 \textrm{ord}=\textrm{'fro'}\end{array}\right. where \sigma_{i} are the singular 
    169 values of \mathbf{A}.  
     182\textrm{ord}=\textrm{'fro'}\end{array}\right</ipython-equation> where <ipython-inlineequation>\sigma_{i}</ipython-inlineequation> are the singular 
     183values of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>.  
    170184</para> 
    171185</section> 
    …… 
    175189<para>Linear least-squares problems occur in many branches of applied 
    176190mathematics. In this problem a set of linear scaling coefficients is sought that 
    177 allow a model to fit data. In particular it is assumed that data y_{i} is 
    178 related to data \mathbf{x}_{i} through a set of coefficients c_{j} and model 
    179 functions f_{j}\left(\mathbf{x}_{i}\right) via the model 
    180 y_{i}=\sum_{j}c_{j}f_{j}\left(\mathbf{x}_{i}\right)+\epsilon_{i} where 
    181 \epsilon_{i} represents uncertainty in the data. The strategy of least squares 
    182 is to pick the coefficients c_{j} to minimize 
    183 J\left(\mathbf{c}\right)=\sum_{i}\left|y_{i}-\sum_{j}c_{j}f_{j}\left(x_{i}\right)\right|^{2}.  
    184 </para> 
    185 <para>Theoretically, a global minimum will occur when \frac{\partial J}{\partial 
    186 c_{n}^{*}}=0=\sum_{i}\left(y_{i}-\sum_{j}c_{j}f_{j}\left(x_{i}\right)\right)\left(-f_{n}^{*}\left(x_{i}\right)\right)or 
    187 \sum_{j}c_{j}\sum_{i}f_{j}\left(x_{i}\right)f_{n}^{*}\left(x_{i}\right) where 
    188 \left\{ \mathbf{A}\right\} _{ij}=f_{j}\left(x_{i}\right). When \mathbf{A^{H}A} 
     191allow a model to fit data. In particular it is assumed that data <ipython-inlineequation>y_{i}</ipython-inlineequation> is 
     192related to data <ipython-inlineequation>\mathbf{x}_{i}</ipython-inlineequation> through a set of coefficients <ipython-inlineequation>c_{j}</ipython-inlineequation> and model 
     193functions <ipython-inlineequation>f_{j}\left(\mathbf{x}_{i}\right)</ipython-inlineequation> via the model 
     194<ipython-equation>y_{i}=\sum_{j}c_{j}f_{j}\left(\mathbf{x}_{i}\right)+\epsilon_{i}</ipython-equation> where 
     195<ipython-inlineequation>\epsilon_{i}</ipython-inlineequation> represents uncertainty in the data. The strategy of least squares 
     196is to pick the coefficients <ipython-inlineequation>c_{j}</ipython-inlineequation> to minimize 
     197<ipython-equation> 
     198J\left(\mathbf{c}\right)=\sum_{i}\left|y_{i}-\sum_{j}c_{j}f_{j}\left(x_{i}\right)\right|^{2}.\] 
     199</ipython-equation></para> 
     200<para>Theoretically, a global minimum will occur when  
     201<ipython-equation>\frac{\partial J}{\partial 
     202c_{n}^{*}}=0=\sum_{i}\left(y_{i}-\sum_{j}c_{j}f_{j}\left(x_{i}\right)\right)\left(-f_{n}^{*}\left(x_{i}\right)\right)</ipython-equation> or 
     203<ipython-equation verb="1">\begin{eqnarray*} 
     204\sum_{j}c_{j}\sum_{i}f_{j}\left(x_{i}\right)f_{n}^{*}\left(x_{i}\right) &amp; = &amp; \sum_{i}y_{i}f_{n}^{*} 
     205\left(x_{i}\right)\\ 
     206\mathbf{A}^{H}\mathbf{Ac} &amp; = &amp; \mathbf{A}^{H}\mathbf{y}\end{eqnarray*} 
     207</ipython-equation> where 
     208<ipython-equation>\left\{ \mathbf{A}\right\} _{ij}=f_{j}\left(x_{i}\right)</ipython-equation>. When <ipython-inlineequation>\mathbf{A^{H}A}</ipython-inlineequation>  
    189209is invertible, then 
    190 \mathbf{c}=\left(\mathbf{A}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H}\mathbf{y}=\mathbf{A}^{\dagger}\mathbf{y} 
    191 where \mathbf{A}^{\dagger} is called the pseudo-inverse of \mathbf{A}. Notice 
    192 that using this definition of \mathbf{A} the model can be written 
    193 \mathbf{y}=\mathbf{Ac}+\boldsymbol{\epsilon}.The command linalg.lstsq will solve 
    194 the linear least squares problem for \mathbf{c} given \mathbf{A} and \mathbf{y}. 
     210<ipython-equation>\mathbf{c}=\left(\mathbf{A}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H}\mathbf{y}=\mathbf{A}^{\dagger}\mathbf{y}</ipython-equation>  
     211where <ipython-inlineequation>\mathbf{A}^{\dagger}</ipython-inlineequation> is called the pseudo-inverse of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>. Notice 
     212that using this definition of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> the model can be written 
     213<ipython-equation>\mathbf{y}=\mathbf{Ac}+\boldsymbol{\epsilon}.</ipython-equation> The command linalg.lstsq will solve 
     214the linear least squares problem for <ipython-inlineequation>\mathbf{c}</ipython-inlineequation> given <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> and <ipython-inlineequation>\mathbf{y}</ipython-inlineequation>. 
    195215In addition linalg.pinv or linalg.pinv2 (uses a different method based on 
    196 singular value decomposition) will find \mathbf{A}^{\dagger} given \mathbf{A}.  
     216singular value decomposition) will find <ipython-inlineequation>\mathbf{A}^{\dagger}</ipython-inlineequation> given <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>.  
    197217</para> 
    198218<para>The following example and figure demonstrate the use of linalg.lstsq and 
    199219linalg.pinv for solving a data-fitting problem. The data shown below were 
    200 generated using the model: y_{i}=c_{1}e^{-x_{i}}+c_{2}x_{i} where x_{i}=0.1i for 
    201 i=1\ldots10, c_{1}=5, and c_{2}=4. Noise is added to y_{i} and the coefficients 
    202 c_{1} and c_{2} are estimated using linear least squares.  
    203 </para> 
    204 <para> XXX: ipython 
     220generated using the model: <ipython-equation>y_{i}=c_{1}e^{-x_{i}}+c_{2}x_{i}</ipython-equation> where <ipython-inlineequation>x_{i}=0.1i</ipython-inlineequation> for 
     221<ipython-inlineequation>i=1\ldots10</ipython-inlineequation>, <ipython-inlineequation>c_{1}=5</ipython-inlineequation>, and <ipython-inlineequation>c_{2}=4</ipython-inlineequation>. Noise is added to <ipython-inlineequation>y_{i}</ipython-inlineequation> and the coefficients 
     222<ipython-inlineequation>c_{1}</ipython-inlineequation> and <ipython-inlineequation>c_{2}</ipython-inlineequation> are estimated using linear least squares.  
     223</para> 
     224<programlisting> XXX: ipython 
    205225c1,c2= 5.0,2.0 
    206226i = r_[1:11] 
    …… 
    220240xplt.title('Data fitting with linalg.lstsq') 
    221241xplt.eps('lstsq_fit') 
    222 </para> 
     242</programlisting> 
    223243<para> 
    224244[Graphics file: lstsq_fit.eps] 
    …… 
    231251or linalg.pinv2. These two commands differ in how they compute the generalized 
    232252inverse. The first uses the linalg.lstsq algorithm while the second uses 
    233 singular value decomposition. Let \mathbf{A} be an M\times N matrix, then if M&gt;N 
     253singular value decomposition. Let <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> be an <ipython-inlineequation>M\times N</ipython-inlineequation> matrix, then if <ipython-inlineequation>M&gt;N</ipython-inlineequation>  
    234254the generalized inverse is 
    235 \mathbf{A}^{\dagger}=\left(\mathbf{A}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H} 
    236 while if M&lt;N matrix the generalized inverse is 
    237 \mathbf{A}^{\#}=\mathbf{A}^{H}\left(\mathbf{A}\mathbf{A}^{H}\right)^{-1}. In 
    238 both cases for M=N, then \mathbf{A}^{\dagger}=\mathbf{A}^{\#}=\mathbf{A}^{-1}as 
    239 long as \mathbf{A} is invertible.  
     255<ipython-equation>\mathbf{A}^{\dagger}=\left(\mathbf{A}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H}</ipython-equation>  
     256while if <ipython-inlineequation>M&lt;N</ipython-inlineequation> matrix the generalized inverse is 
     257<ipython-equation>\mathbf{A}^{\#}=\mathbf{A}^{H}\left(\mathbf{A}\mathbf{A}^{H}\right)^{-1}.</ipython-equation> In 
     258both cases for <ipython-inlineequation>M=N</ipython-inlineequation>, then <ipython-equation>\mathbf{A}^{\dagger}=\mathbf{A}^{\#}=\mathbf{A}^{-1}</ipython-equation> as 
     259long as <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> is invertible.  
    240260</para> 
    241261</section> 
    …… 
    252272<para>The eigenvalue-eigenvector problem is one of the most commonly 
    253273employed linear algebra operations. In one popular form, the 
    254 eigenvalue-eigenvector problem is to find for some square matrix \mathbf{A} 
    255 scalars \lambda  and corresponding vectors \mathbf{v} such that 
    256 \mathbf{Av}=\lambda\mathbf{v}. For an N\times N matrix, there are N (not 
     274eigenvalue-eigenvector problem is to find for some square matrix <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>  
     275scalars <ipython-inlineequation>\lambda</ipython-inlineequation>  and corresponding vectors <ipython-inlineequation>\mathbf{v}</ipython-inlineequation> such that 
     276<ipython-equation>\mathbf{Av}=\lambda\mathbf{v}.</ipython-equation> For an <ipython-inlineequation>N\times N</ipython-inlineequation> matrix, there are N (not 
    257277necessarily distinct) eigenvalues --- roots of the (characteristic) polynomial 
    258 \left|\mathbf{A}-\lambda\mathbf{I}\right|=0. 
    259 </para> 
    260 <para>The eigenvectors, \mathbf{v}, are also sometimes called right eigenvectors to 
     278<ipython-equation>\left|\mathbf{A}-\lambda\mathbf{I}\right|=0.</ipython-equation>  
     279</para> 
     280<para>The eigenvectors, <ipython-inlineequation>\mathbf{v}</ipython-inlineequation>, are also sometimes called right eigenvectors to 
    261281distinguish them from another set of left eigenvectors that satisfy 
    262 \mathbf{v}_{L}^{H}\mathbf{A}=\lambda\mathbf{v}_{L}^{H}or 
    263 \mathbf{A}^{H}\mathbf{v}_{L}=\lambda^{*}\mathbf{v}_{L}. With it's default 
    264 optional arguments, the command linalg.eig returns \lambda  and \mathbf{v}. 
    265 However, it can also return \mathbf{v}_{L} and just \lambda  by itself 
    266 (linalg.eigvals returns just \lambda  as well).  
     282<ipython-equation>\mathbf{v}_{L}^{H}\mathbf{A}=\lambda\mathbf{v}_{L}^{H}</ipython-equation> or 
     283<ipython-equation>\mathbf{A}^{H}\mathbf{v}_{L}=\lambda^{*}\mathbf{v}_{L}.</ipython-equation> With it's default 
     284optional arguments, the command linalg.eig returns <ipython-inlineequation>\lambda</ipython-inlineequation>  and <ipython-inlineequation>\mathbf{v}</ipython-inlineequation>. 
     285However, it can also return <ipython-inlineequation>\mathbf{v}_{L}</ipython-inlineequation> and just <ipython-inlineequation>\lambda</ipython-inlineequation>  by itself 
     286(linalg.eigvals returns just <ipython-inlineequation>\lambda</ipython-inlineequation>  as well).  
    267287</para> 
    268288<para>In addtion, linalg.eig can also solve the more general eigenvalue problem 
    269 \mathbf{Av}for square matrices \mathbf{A} and \mathbf{B}. The standard 
     289<ipython-equation verb="1">\begin{eqnarray*} 
     290\mathbf{Av} &amp; = &amp; \lambda\mathbf{Bv}\\ 
     291\mathbf{A}^{H}\mathbf{v}_{L} &amp; = &amp; \lambda^{*}\mathbf{B}^{H}\mathbf{v}_{L}\end{eqnarray*}</ipython-equation> 
     292for square matrices <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> and <ipython-inlineequation>\mathbf{B}</ipython-inlineequation>. The standard 
    270293eigenvalue problem is an example of the general eigenvalue problem for 
    271 \mathbf{B}=\mathbf{I}. When a generalized eigenvalue problem can be solved, then 
    272 it provides a decomposition of \mathbf{A} as 
    273 \mathbf{A}=\mathbf{BV}\boldsymbol{\Lambda}\mathbf{V}^{-1} where \mathbf{V} is 
    274 the collection of eigenvectors into columns and \boldsymbol{\Lambda} is a 
     294<ipython-inlineequation>\mathbf{B}=\mathbf{I}</ipython-inlineequation>. When a generalized eigenvalue problem can be solved, then 
     295it provides a decomposition of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> as 
     296<ipython-equation>\mathbf{A}=\mathbf{BV}\boldsymbol{\Lambda}\mathbf{V}^{-1}</ipython-equation> where <ipython-inlineequation>\mathbf{V}</ipython-inlineequation> is 
     297the collection of eigenvectors into columns and <ipython-inlineequation>\boldsymbol{\Lambda}</ipython-inlineequation> is a 
    275298diagonal matrix of eigenvalues.  
    276299</para> 
    277300<para>By definition, eigenvectors are only defined up to a constant scale factor. In 
    278 SciPy, the scaling factor for the eigenvectors is chosen so that \left\Vert 
    279 \mathbf{v}\right\Vert ^{2}=\sum_{i}v_{i}^{2}=1.  
     301SciPy, the scaling factor for the eigenvectors is chosen so that <ipython-inlineequation>\left\Vert 
     302\mathbf{v}\right\Vert ^{2}=\sum_{i}v_{i}^{2}=1</ipython-inlineequation>.  
    280303</para> 
    281304<para>As an example, consider finding the eigenvalues and eigenvectors of the matrix 
    282 \mathbf{A}=\left[\begin{array}{ccc} 
     305<ipython-equation>\mathbf{A}=\left[\begin{array}{ccc} 
    2833061 &amp; 5 &amp; 2\\ 
    2843072 &amp; 4 &amp; 1\\ 
    285 3 &amp; 6 &amp; 2\end{array}\right]. The characteristic polynomial is 
    286 \left|\mathbf{A}-\lambda\mathbf{I}\right| The roots of this polynomial are the 
    287 eigenvalues of \mathbf{A}: \lambda_{1} The eigenvectors corresponding to each 
     3083 &amp; 6 &amp; 2\end{array}\right].</ipython-equation> The characteristic polynomial is 
     309<ipython-equation verb="1"> 
     310\begin{eqnarray*} 
     311\left|\mathbf{A}-\lambda\mathbf{I}\right| &amp; = &amp; \left(1-\lambda\right)\left[\left(4-\lambda\right)\left(2-\lambda\right)-6\right]-\\ 
     312 &amp;  &amp; 5\left[2\left(2-\lambda\right)-3\right]+2\left[12-3\left(4-\lambda\right)\right]\\ 
     313 &amp; = &amp; -\lambda^{3}+7\lambda^{2}+8\lambda-3.\end{eqnarray*} 
     314</ipython-equation> 
     315 The roots of this polynomial are the 
     316eigenvalues of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>:  
     317<ipython-equation verb="1"> 
     318\begin{eqnarray*} 
     319\lambda_{1} &amp; = &amp; 7.9579\\ 
     320\lambda_{2} &amp; = &amp; -1.2577\\ 
     321\lambda_{3} &amp; = &amp; 0.2997.\end{eqnarray*} 
     322</ipython-equation> 
     323 The eigenvectors corresponding to each 
    288324eigenvalue can be found using the original equation. The eigenvectors associated 
    289325with these eigenvalues can then be found.  
    290326</para> 
    291 <para>XXX: ipython 
     327<programlisting>XXX: ipython 
    292328&gt;&gt;&gt; A = mat('[1 5 2; 2 4 1; 3 6 2]') 
    293329&gt;&gt;&gt; la,v = linalg.eig(A) 
    …… 
    308344&gt;&gt;&gt; print max(ravel(abs(A*v1-l1*v1))) 
    3093454.4408920985e-16 
    310 </para> 
     346</programlisting> 
    311347</section> 
    312348 
    …… 
    314350<title>Singular value decomposition</title> 
    315351<para>Singular Value Decompostion (SVD) can be thought of as an extension of the 
    316 eigenvalue problem to matrices that are not square. Let \mathbf{A} be an M\times 
    317 N matrix with M and N arbitrary. The matrices \mathbf{A}^{H}\mathbf{A} and 
    318 \mathbf{A}\mathbf{A}^{H} are square hermitian matrices<footnote><para>A hermition matrix 
    319 \mathbf{D} satisfies \mathbf{D}^{H}=\mathbf{D}.</para></footnote> of size N\times N and M\times M 
     352eigenvalue problem to matrices that are not square. Let <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> be an  
     353<ipython-inlineequation>M\times N</ipython-inlineequation> matrix with <ipython-inlineequation>M</ipython-inlineequation> and <ipython-inlineequation>N</ipython-inlineequation> arbitrary. The matrices <ipython-inlineequation>\mathbf{A}^{H}\mathbf{A}</ipython-inlineequation> and 
     354<ipython-inlineequation>\mathbf{A}\mathbf{A}^{H}</ipython-inlineequation> are square hermitian matrices<footnote><para>A hermition matrix 
     355<ipython-inlineequation>\mathbf{D}</ipython-inlineequation> satisfies <ipython-inlineequation>\mathbf{D}^{H}=\mathbf{D}</ipython-inlineequation>.</para></footnote> of size <ipython-inlineequation>N\times N</ipython-inlineequation> and <ipython-inlineequation>M\times M</ipython-inlineequation>  
    320356respectively. It is known that the eigenvalues of square hermitian matrices are 
    321 real and non-negative. In addtion, there are at most \min\left(M,N\right) 
    322 identical non-zero eigenvalues of \mathbf{A}^{H}\mathbf{A} and 
    323 \mathbf{A}\mathbf{A}^{H}. Define these positive eigenvalues as \sigma_{i}^{2}. 
    324 The square-root of these are called singular values of \mathbf{A}. The 
    325 eigenvectors of \mathbf{A}^{H}\mathbf{A} are collected by columns into an 
    326 N\times N unitary<footnote><para>A unitary matrix \mathbf{D} satisfies 
    327 \mathbf{D}^{H}\mathbf{D}=\mathbf{I}=\mathbf{D}\mathbf{D}^{H} so that 
    328 \mathbf{D}^{-1}=\mathbf{D}^{H}.</para></footnote> matrix \mathbf{V} while the eigenvectors of 
    329 \mathbf{A}\mathbf{A}^{H} are collected by columns in the unitary matrix 
    330 \mathbf{U}, the singular values are collected in an M\times N zero matrix 
    331 \mathbf{\boldsymbol{\Sigma}} with main diagonal entries set to the singular 
    332 values. Then \mathbf{A=U}\boldsymbol{\Sigma}\mathbf{V}^{H} is the singular-value 
    333 decomposition of \mathbf{A}. Every matrix has a singular value decomposition. 
    334 Sometimes, the singular values are called the spectrum of \mathbf{A}. The 
    335 command linalg.svd will return \mathbf{U}, \mathbf{V}^{H}, and \sigma_{i} as an 
    336 array of the singular values. To obtain the matrix \mathbf{\Sigma} use 
     357real and non-negative. In addtion, there are at most <ipython-inlineequation>\min\left(M,N\right)</ipython-inlineequation>  
     358identical non-zero eigenvalues of <ipython-inlineequation>\mathbf{A}^{H}\mathbf{A}</ipython-inlineequation> and 
     359<ipython-inlineequation>\mathbf{A}\mathbf{A}^{H}</ipython-inlineequation>. Define these positive eigenvalues as <ipython-inlineequation>\sigma_{i}^{2}</ipython-inlineequation>. 
     360The square-root of these are called singular values of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>. The 
     361eigenvectors of <ipython-inlineequation>\mathbf{A}^{H}\mathbf{A}</ipython-inlineequation> are collected by columns into an 
     362<ipython-inlineequation>N\times N</ipython-inlineequation> unitary<footnote><para>A unitary matrix <ipython-inlineequation>\mathbf{D}</ipython-inlineequation> satisfies 
     363<ipython-inlineequation>\mathbf{D}^{H}\mathbf{D}=\mathbf{I}=\mathbf{D}\mathbf{D}^{H}</ipython-inlineequation> so that 
     364<ipython-inlineequation>\mathbf{D}^{-1}=\mathbf{D}^{H}</ipython-inlineequation>.</para></footnote> matrix <ipython-inlineequation>\mathbf{V}</ipython-inlineequation> while the eigenvectors of 
     365<ipython-inlineequation>\mathbf{A}\mathbf{A}^{H}</ipython-inlineequation> are collected by columns in the unitary matrix 
     366<ipython-inlineequation>\mathbf{U}</ipython-inlineequation>, the singular values are collected in an <ipython-inlineequation>M\times N</ipython-inlineequation> zero matrix 
     367<ipython-inlineequation>\mathbf{\boldsymbol{\Sigma}}</ipython-inlineequation> with main diagonal entries set to the singular 
     368values. Then <ipython-inlineequation>\mathbf{A=U}\boldsymbol{\Sigma}\mathbf{V}^{H}</ipython-inlineequation> is the singular-value 
     369decomposition of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>. Every matrix has a singular value decomposition. 
     370Sometimes, the singular values are called the spectrum of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>. The 
     371command linalg.svd will return <ipython-inlineequation>\mathbf{U}</ipython-inlineequation>, <ipython-inlineequation>\mathbf{V}^{H}</ipython-inlineequation>, and <ipython-inlineequation>\sigma_{i}</ipython-inlineequation> as an 
     372array of the singular values. To obtain the matrix <ipython-inlineequation>\mathbf{\Sigma}</ipython-inlineequation> use 
    337373linalg.diagsvd. The following example illustrates the use of linalg.svd.  
    338374</para> 
    339 <para>XXX: ipython 
     375<programlisting>XXX: ipython 
    340376&gt;&gt;&gt; A = mat('[1 3 2; 1 2 3]') 
    341377&gt;&gt;&gt; M,N = A.shape 
    …… 
    360396Matrix([[ 1.,  3.,  2.], 
    361397       [ 1.,  2.,  3.]]) 
    362 </para> 
     398</programlisting> 
    363399</section> 
    364400 
    365401<section> 
    366402<title>LU decomposition</title> 
    367 <para>The LU decompostion finds a representation for the M\times N matrix 
    368 \mathbf{A} as \mathbf{A}=\mathbf{PLU} where \mathbf{P} is an M\times M 
     403<para>The LU decompostion finds a representation for the <ipython-inlineequation>M\times N</ipython-inlineequation> matrix 
     404<ipython-inlineequation>\mathbf{A}</ipython-inlineequation> as <ipython-inlineequation>\mathbf{A}=\mathbf{PLU}</ipython-inlineequation> where <ipython-inlineequation>\mathbf{P}</ipython-inlineequation> is an <ipython-inlineequation>M\times M</ipython-inlineequation>  
    369405permutation matrix (a permutation of the rows of the identity matrix), 
    370 \mathbf{L} is in M\times K lower triangular or trapezoidal matrix 
    371 (K=\min\left(M,N\right)) with unit-diagonal, and \mathbf{U} is an upper 
     406<ipython-inlineequation>\mathbf{L}</ipython-inlineequation> is in <ipython-inlineequation>M\times K</ipython-inlineequation> lower triangular or trapezoidal matrix 
     407(<ipython-inlineequation>K=\min\left(M,N\right)</ipython-inlineequation>) with unit-diagonal, and <ipython-inlineequation>\mathbf{U}</ipython-inlineequation> is an upper 
    372408triangular or trapezoidal matrix. The SciPy command for this decomposition is 
    373409linalg.lu.  
    …… 
    376412Such a decomposition is often useful for solving many simultaneous equations 
    377413where the left-hand-side does not change but the right hand side does. For 
    378 example, suppose we are going to solve \mathbf{A}\mathbf{x}_{i}=\mathbf{b}_{i} 
    379 for many different \mathbf{b}_{i}. The LU decomposition allows this to be 
    380 written as \mathbf{PLUx}_{i}=\mathbf{b}_{i}. Because \mathbf{L} is 
    381 lower-triangular, the equation can be solved for \mathbf{U}\mathbf{x}_{i} and 
    382 finally \mathbf{x}_{i} very rapidly using forward- and back-substitution. An 
    383 initial time spent factoring \mathbf{A} allows for very rapid solution of 
     414example, suppose we are going to solve <ipython-inlineequation>\mathbf{A}\mathbf{x}_{i}=\mathbf{b}_{i}</ipython-inlineequation>  
     415for many different <ipython-inlineequation>\mathbf{b}_{i}</ipython-inlineequation>. The LU decomposition allows this to be 
     416written as <ipython-inlineequation>\mathbf{PLUx}_{i}=\mathbf{b}_{i}</ipython-inlineequation>. Because <ipython-inlineequation>\mathbf{L}</ipython-inlineequation> is 
     417lower-triangular, the equation can be solved for <ipython-inlineequation>\mathbf{U}\mathbf{x}_{i}</ipython-inlineequation> and 
     418finally <ipython-inlineequation>\mathbf{x}_{i}</ipython-inlineequation> very rapidly using forward- and back-substitution. An 
     419initial time spent factoring <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> allows for very rapid solution of 
    384420similar systems of equations in the future. If the intent for performing LU 
    385421decomposition is for solving linear systems then the command linalg.lu_factor 
    …… 
    393429<para>Cholesky decomposition is a special case of LU decomposition 
    394430applicable to Hermitian positive definite matrices. When 
    395 \mathbf{A}=\mathbf{A}^{H} and \mathbf{x}^{H}\mathbf{Ax}\geq0 for all \mathbf{x}, 
    396 then decompositions of \mathbf{A} can be found so that \mathbf{A}where 
    397 \mathbf{L} is lower-triangular and \mathbf{U} is upper triangular. Notice that 
    398 \mathbf{L}=\mathbf{U}^{H}. The command linagl.cholesky computes the cholesky 
     431<ipython-inlineequation>\mathbf{A}=\mathbf{A}^{H}</ipython-inlineequation> and <ipython-inlineequation>\mathbf{x}^{H}\mathbf{Ax}\geq0</ipython-inlineequation> for all <ipython-inlineequation>\mathbf{x}</ipython-inlineequation>, 
     432then decompositions of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> can be found so that 
     433<ipython-equation verb="1"> 
     434\begin{eqnarray*} 
     435\mathbf{A} &amp; = &amp; \mathbf{U}^{H}\mathbf{U}\\ 
     436\mathbf{A} &amp; = &amp; \mathbf{L}\mathbf{L}^{H}\end{eqnarray*} 
     437</ipython-equation> 
     438where 
     439<ipython-inlineequation>\mathbf{L}</ipython-inlineequation> is lower-triangular and <ipython-inlineequation>\mathbf{U}</ipython-inlineequation> is upper triangular. Notice that 
     440<ipython-inlineequation>\mathbf{L}=\mathbf{U}^{H}</ipython-inlineequation>. The command linagl.cholesky computes the cholesky 
    399441factorization. For using cholesky factorization to solve systems of equations 
    400442there are also linalg.cho_factor and linalg.cho_solve routines that work 
    …… 
    406448<title>QR decomposition</title> 
    407449<para>The QR decomposition (sometimes called a polar decomposition) works for 
    408 any M\times N array and finds an M\times M unitary matrix \mathbf{Q} and an 
    409 M\times N upper-trapezoidal matrix \mathbf{R} such that \mathbf{A=QR}. Notice 
    410 that if the SVD of \mathbf{A} is known then the QR decomposition can be found 
    411 \mathbf{A}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{H}=\mathbf{QR} implies that 
    412 \mathbf{Q}=\mathbf{U} and \mathbf{R}=\boldsymbol{\Sigma}\mathbf{V}^{H}. Note, 
     450any <ipython-inlineequation>M\times N</ipython-inlineequation> array and finds an <ipython-inlineequation>M\times M</ipython-inlineequation> unitary matrix <ipython-inlineequation>\mathbf{Q}</ipython-inlineequation> and an 
     451<ipython-inlineequation>M\times N</ipython-inlineequation> upper-trapezoidal matrix <ipython-inlineequation>\mathbf{R}</ipython-inlineequation> such that <ipython-inlineequation>\mathbf{A=QR}</ipython-inlineequation>. Notice 
     452that if the SVD of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> is known then the QR decomposition can be found 
     453<ipython-equation>\mathbf{A}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{H}=\mathbf{QR}</ipython-equation> implies that 
     454<ipython-inlineequation>\mathbf{Q}=\mathbf{U}</ipython-inlineequation> and <ipython-inlineequation>\mathbf{R}=\boldsymbol{\Sigma}\mathbf{V}^{H}</ipython-inlineequation>. Note, 
    413455however, that in SciPy independent algorithms are used to find QR and SVD 
    414456decompositions. The command for QR decomposition is linalg.qr.  
    …… 
    418460<section> 
    419461<title>Schur decomposition</title> 
    420 <para>For a square N\times N matrix, \mathbf{A}, the Schur decomposition finds 
    421 (not-necessarily unique) matrices \mathbf{T} and \mathbf{Z} such that 
    422 \mathbf{A}=\mathbf{ZT}\mathbf{Z}^{H} where \mathbf{Z} is a unitary matrix and 
    423 \mathbf{T} is either upper-triangular or quasi-upper triangular depending on 
     462<para>For a square <ipython-inlineequation>N\times N</ipython-inlineequation> matrix, <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>, the Schur decomposition finds 
     463(not-necessarily unique) matrices <ipython-inlineequation>\mathbf{T}</ipython-inlineequation> and <ipython-inlineequation>\mathbf{Z}</ipython-inlineequation> such that 
     464<ipython-equation>\mathbf{A}=\mathbf{ZT}\mathbf{Z}^{H}</ipython-equation> where <ipython-inlineequation>\mathbf{Z}</ipython-inlineequation> is a unitary matrix and 
     465<ipython-inlineequation>\mathbf{T}</ipython-inlineequation> is either upper-triangular or quasi-upper triangular depending on 
    424466whether or not a real schur form or complex schur form is requested. For a real 
    425 schur form both \mathbf{T} and \mathbf{Z} are real-valued when \mathbf{A} is 
    426 real-valued. When \mathbf{A} is a real-valued matrix the real schur form is only 
    427 quasi-upper triangular because 2\times2 blocks extrude from the main diagonal 
     467schur form both <ipython-inlineequation>\mathbf{T}</ipython-inlineequation> and <ipython-inlineequation>\mathbf{Z}</ipython-inlineequation> are real-valued when <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> is 
     468real-valued. When <ipython-inlineequation>\mathbf{A}</ipython-inlineequation> is a real-valued matrix the real schur form is only 
     469quasi-upper triangular because <ipython-inlineequation>2\times2</ipython-inlineequation> blocks extrude from the main diagonal 
    428470corresponding to any complex-valued eigenvalues. The command linalg.schur finds 
    429 the Schur decomposition while the command linalg.rsf2csf converts \mathbf{T} and 
    430 \mathbf{Z} from a real Schur form to a complex Schur form. The Schur form is 
     471the Schur decomposition while the command linalg.rsf2csf converts <ipython-inlineequation>\mathbf{T}</ipython-inlineequation> and 
     472<ipython-inlineequation>\mathbf{Z}</ipython-inlineequation> from a real Schur form to a complex Schur form. The Schur form is 
    431473especially useful in calculating functions of matrices.  
    432474</para> 
    433475<para>The following example illustrates the schur decomposition: 
    434476</para> 
    435 <para>XXX: ipython 
     477<programlisting>XXX: ipython 
    436478&gt;&gt;&gt; A = mat('[1 3 2; 1 4 5; 2 3 6]') 
    437479&gt;&gt;&gt; T,Z = linalg.schur(A) 
    …… 
    467509       [ 0.,  0.,  0.], 
    468510       [ 0.,  0.,  0.]]) 
    469 </para> 
     511</programlisting> 
    470512</section> 
    471513</section> 
    …… 
    473515<section> 
    474516<title>Matrix Functions</title> 
    475 <para>Consider the function f\left(x\right) with Taylor series expansion 
    476 f\left(x\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}x^{k}. 
     517<para>Consider the function <ipython-inlineequation>f\left(x\right)</ipython-inlineequation> with Taylor series expansion 
     518<ipython-equation>f\left(x\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}x^{k}.</ipython-equation>  
    477519A matrix function can be defined using this Taylor series for the square matrix 
    478 \mathbf{A} as 
    479 f\left(\mathbf{A}\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}\mathbf{A}^{k}. 
     520<ipython-inlineequation>\mathbf{A}</ipython-inlineequation> as 
     521<ipython-equation>f\left(\mathbf{A}\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}\mathbf{A}^{k}.</ipython-equation>  
    480522While, this serves as a useful representation of a matrix function, it is rarely 
    481523the best way to calculate a matrix function.  
    …… 
    486528<para>The matrix exponential is one of the more common matrix functions. It can 
    487529be defined for square matrices as 
    488 e^{\mathbf{A}}=\sum_{k=0}^{\infty}\frac{1}{k!}\mathbf{A}^{k}. The command 
     530<ipython-equation>e^{\mathbf{A}}=\sum_{k=0}^{\infty}\frac{1}{k!}\mathbf{A}^{k}.</ipython-equation> The command 
    489531linalg.expm3 uses this Taylor series definition to compute the matrix 
    490532exponential. Due to poor convergence properties it is not often used.  
    491533</para> 
    492534<para>Another method to compute the matrix exponential is to find an eigenvalue 
    493 decomposition of \mathbf{A}: 
    494 \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} and note that 
    495 e^{\mathbf{A}}=\mathbf{V}e^{\boldsymbol{\Lambda}}\mathbf{V}^{-1} where the 
    496 matrix exponential of the diagonal matrix \boldsymbol{\Lambda} is just the 
     535decomposition of <ipython-inlineequation>\mathbf{A}</ipython-inlineequation>: 
     536<ipython-equation>\mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}</ipython-equation> and note that 
     537<ipython-equation>e^{\mathbf{A}}=\mathbf{V}e^{\boldsymbol{\Lambda}}\mathbf{V}^{-1}</ipython-equation> where the 
     538matrix exponential of the diagonal matrix <ipython-inlineequation>\boldsymbol{\Lambda}</ipython-inlineequation> is just the 
    497539exponential of its elements. This method is implemented in linalg.expm2.  
    498540</para> 
    499541<para>The preferred method for implementing the matrix exponential is to use scaling 
    500 and a Padé approximation for e^{x}. This algorithm is implemented as 
     542and a Padé approximation for <ipython-inlineequation&g