# 线性代数

`numpy` 和 `scipy` 中，负责进行线性代数部分计算的模块叫做 `linalg`。

In [1]:

```py
import numpy as np
import numpy.linalg
import scipy as sp
import scipy.linalg
import matplotlib.pyplot as plt
from scipy import linalg

%matplotlib inline

```

## numpy.linalg VS scipy.linalg

一方面`scipy.linalg` 包含 `numpy.linalg` 中的所有函数，同时还包含了很多 `numpy.linalg` 中没有的函数。

另一方面，`scipy.linalg` 能够保证这些函数使用 BLAS/LAPACK 加速，而 `numpy.linalg` 中这些加速是可选的。

因此，在使用时，我们一般使用 `scipy.linalg` 而不是 `numpy.linalg`。

我们可以简单看看两个模块的差异：

In [2]:

```py
print "number of items in numpy.linalg:", len(dir(numpy.linalg))
print "number of items in scipy.linalg:", len(dir(scipy.linalg))

```

```py
number of items in numpy.linalg: 36
number of items in scipy.linalg: 115

```

## numpy.matrix VS 2D numpy.ndarray

线性代数的基本操作对象是矩阵，而矩阵的表示方法主要有两种：`numpy.matrix` 和 2D `numpy.ndarray`。

### numpy.matrix

`numpy.matrix` 是一个矩阵类，提供了一些方便的矩阵操作：

*   支持类似 `MATLAB` 创建矩阵的语法
*   矩阵乘法默认用 `*` 号
*   `.I` 表示逆，`.T` 表示转置

可以用 `mat` 或者 `matrix` 来产生矩阵：

In [3]:

```py
A = np.mat("[1, 2; 3, 4]")
print repr(A)

A = np.matrix("[1, 2; 3, 4]")
print repr(A)

```

```py
matrix([[1, 2],
        [3, 4]])
matrix([[1, 2],
        [3, 4]])

```

转置和逆：

In [4]:

```py
print repr(A.I)
print repr(A.T)

```

```py
matrix([[-2\. ,  1\. ],
        [ 1.5, -0.5]])
matrix([[1, 3],
        [2, 4]])

```

矩阵乘法：

In [5]:

```py
b = np.mat('[5; 6]')
print repr(A * b)

```

```py
matrix([[17],
        [39]])

```

### 2 维 numpy.ndarray

虽然 `numpy.matrix` 有着上面的好处，但是一般不建议使用，而是用 2 维 `numpy.ndarray` 对象替代，这样可以避免一些不必要的困惑。

我们可以使用 `array` 复现上面的操作：

In [6]:

```py
A = np.array([[1,2], [3,4]])
print repr(A)

```

```py
array([[1, 2],
       [3, 4]])

```

逆和转置：

In [7]:

```py
print repr(linalg.inv(A))
print repr(A.T)

```

```py
array([[-2\. ,  1\. ],
       [ 1.5, -0.5]])
array([[1, 3],
       [2, 4]])

```

矩阵乘法：

In [8]:

```py
b = np.array([5, 6])

print repr(A.dot(b))

```

```py
array([17, 39])

```

普通乘法：

In [9]:

```py
print repr(A * b)

```

```py
array([[ 5, 12],
       [15, 24]])

```

`scipy.linalg` 的操作可以作用到两种类型的对象上，没有区别。

## 基本操作

### 求逆

矩阵 $\mathbf{A}$ 的逆 $\mathbf{B}$ 满足：$\mathbf{BA}=\mathbf{AB}=I$，记作 $\mathbf{B} = \mathbf{A}^{-1}$。

事实上，我们已经见过求逆的操作，`linalg.inv` 可以求一个可逆矩阵的逆：

In [10]:

```py
A = np.array([[1,2],[3,4]])

print linalg.inv(A)

print A.dot(scipy.linalg.inv(A))

```

```py
[[-2\.   1\. ]
 [ 1.5 -0.5]]
[[  1.00000000e+00   0.00000000e+00]
 [  8.88178420e-16   1.00000000e+00]]

```

### 求解线性方程组

例如，下列方程组 $$ \begin{eqnarray*} x + 3y + 5z & = & 10 \\ 2x + 5y + z & = & 8 \\ 2x + 3y + 8z & = & 3 \end{eqnarray*} $$ 的解为： $$ \begin{split}\left[\begin{array}{c} x\\ y\\ z\end{array}\right]=\left[\begin{array}{ccc} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8\end{array}\right]^{-1}\left[\begin{array}{c} 10\\ 8\\ 3\end{array}\right]=\frac{1}{25}\left[\begin{array}{c} -232\\ 129\\ 19\end{array}\right]=\left[\begin{array}{c} -9.28\\ 5.16\\ 0.76\end{array}\right].\end{split} $$

我们可以使用 `linalg.solve` 求解方程组，也可以先求逆再相乘，两者中 `solve` 比较快。

In [11]:

```py
import time

A = np.array([[1, 3, 5],
              [2, 5, 1],
              [2, 3, 8]])
b = np.array([10, 8, 3])

tic = time.time()

for i in xrange(1000):
    x = linalg.inv(A).dot(b)

print x
print A.dot(x)-b
print "inv and dot: {} s".format(time.time() - tic)

tic = time.time()

for i in xrange(1000):
    x = linalg.solve(A, b)

print x
print A.dot(x)-b
print "solve: {} s".format(time.time() - tic)

```

```py
[-9.28  5.16  0.76]
[  0.00000000e+00  -1.77635684e-15  -8.88178420e-16]
inv and dot: 0.0353579521179 s
[-9.28  5.16  0.76]
[  0.00000000e+00  -1.77635684e-15  -1.77635684e-15]
solve: 0.0284671783447 s

```

### 计算行列式

方阵的行列式为 $$ \left|\mathbf{A}\right|=\sum_{j}\left(-1\right)^{i+j}a_{ij}M_{ij}. $$

其中 $a_{ij}$ 表示 $\mathbf{A}$ 的第 $i$ 行 第 $j$ 列的元素，$M_{ij}$ 表示矩阵 $\mathbf{A}$ 去掉第 $i$ 行 第 $j$ 列的新矩阵的行列式。

例如，矩阵 $$ \begin{split}\mathbf{A=}\left[\begin{array}{ccc} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8\end{array}\right]\end{split} $$ 的行列式是： $$ \begin{eqnarray*} \left|\mathbf{A}\right| & = & 1\left|\begin{array}{cc} 5 & 1\\ 3 & 8\end{array}\right|-3\left|\begin{array}{cc} 2 & 1\\ 2 & 8\end{array}\right|+5\left|\begin{array}{cc} 2 & 5\\ 2 & 3\end{array}\right|\\ & = & 1\left(5\cdot8-3\cdot1\right)-3\left(2\cdot8-2\cdot1\right)+5\left(2\cdot3-2\cdot5\right)=-25.\end{eqnarray*} $$

可以用 `linalg.det` 计算行列式：

In [12]:

```py
A = np.array([[1, 3, 5],
              [2, 5, 1],
              [2, 3, 8]])

print linalg.det(A)

```

```py
-25.0

```

### 计算矩阵或向量的模

矩阵的模定义如下： $$ \begin{split}\left\Vert \mathbf{A}\right\Vert =\left\{ \begin{array}{cc} \max_{i}\sum_{j}\left|a_{ij}\right| & \textrm{ord}=\textrm{inf}\\ \min_{i}\sum_{j}\left|a_{ij}\right| & \textrm{ord}=-\textrm{inf}\\ \max_{j}\sum_{i}\left|a_{ij}\right| & \textrm{ord}=1\\ \min_{j}\sum_{i}\left|a_{ij}\right| & \textrm{ord}=-1\\ \max\sigma_{i} & \textrm{ord}=2\\ \min\sigma_{i} & \textrm{ord}=-2\\ \sqrt{\textrm{trace}\left(\mathbf{A}^{H}\mathbf{A}\right)} & \textrm{ord}=\textrm{'fro'}\end{array}\right.\end{split} $$ 其中，$\sigma_i$ 是矩阵的奇异值。

向量的模定义如下： $$ \begin{split}\left\Vert \mathbf{x}\right\Vert =\left\{ \begin{array}{cc} \max\left|x_{i}\right| & \textrm{ord}=\textrm{inf}\\ \min\left|x_{i}\right| & \textrm{ord}=-\textrm{inf}\\ \left(\sum_{i}\left|x_{i}\right|^{\textrm{ord}}\right)^{1/\textrm{ord}} & \left|\textrm{ord}\right|<\infty.\end{array}\right.\end{split} $$<="" p="">

`linalg.norm` 可以计算向量或者矩阵的模：

In [13]:

```py
A = np.array([[1, 2],
              [3, 4]])

print linalg.norm(A)

print linalg.norm(A,'fro') # frobenius norm 默认值

print linalg.norm(A,1) # L1 norm 最大列和

print linalg.norm(A,-1) # L -1 norm 最小列和

print linalg.norm(A,np.inf) # L inf norm 最大行和

```

```py
5.47722557505
5.47722557505
6
4
7

```

### 最小二乘解和伪逆

#### 问题描述

所谓最小二乘问题的定义如下：

假设 $y_i$ 与 $\mathbf{x_i}$ 的关系可以用一组系数 $c_j$ 和对应的模型函数 $f_j(\mathbf{x_i})$ 的模型表示：

$$ y_{i}=\sum_{j}c_{j}f_{j}\left(\mathbf{x}_{i}\right)+\epsilon_{i} $$

其中 $\epsilon_i$ 表示数据的不确定性。最小二乘就是要优化这样一个关于 $c_j$ 的问题： $$ J\left(\mathbf{c}\right)=\sum_{i}\left|y_{i}-\sum_{j}c_{j}f_{j}\left(x_{i}\right)\right|^{2} $$

其理论解满足： $$ \frac{\partial J}{\partial 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) $$

改写为： $$ \begin{eqnarray*} \sum_{j}c_{j}\sum_{i}f_{j}\left(x_{i}\right)f_{n}^{*}\left(x_{i}\right) & = & \sum_{i}y_{i}f_{n}^{*}\left(x_{i}\right)\\ \mathbf{A}^{H}\mathbf{Ac} & = & \mathbf{A}^{H}\mathbf{y}\end{eqnarray*} $$

其中： $$ \left\{ \mathbf{A}\right\} _{ij}=f_{j}\left(x_{i}\right). $$

当 $\mathbf{A^HA}$ 可逆时，我们有： $$ \mathbf{c}=\left(\mathbf{A}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H}\mathbf{y}=\mathbf{A}^{\dagger}\mathbf{y} $$

矩阵 $\mathbf{A}^{\dagger}$ 叫做 $\mathbf{A}$ 的伪逆。

#### 问题求解

注意到，我们的模型可以写为： $$ \mathbf{y}=\mathbf{Ac}+\boldsymbol{\epsilon}. $$

在给定 $\mathbf{y}$ 和 $\mathbf{A}$ 的情况下，我们可以使用 `linalg.lstsq` 求解 $\mathbf c$。

在给定 $\mathbf{A}$ 的情况下，我们可以使用 `linalg.pinv` 或者 `linalg.pinv2` 求解 $\mathbf{A}^{\dagger}$。

#### 例子

假设我们的数据满足： $$ \begin{align} y_{i} & =c_{1}e^{-x_{i}}+c_{2}x_{i} \\ z_{i} & = y_i + \epsilon_i \end{align} $$

其中 $x_i = \frac{i}{10},\ i = 1,\dots,10$，$c_1 = 5, c_2 = 2$，产生数据

In [14]:

```py
c1, c2 = 5.0, 2.0
i = np.r_[1:11]
xi = 0.1*i
yi = c1*np.exp(-xi) + c2*xi
zi = yi + 0.05 * np.max(yi) * np.random.randn(len(yi))

```

构造矩阵 $\mathbf A$：

In [15]:

```py
A = np.c_[np.exp(-xi)[:, np.newaxis], xi[:, np.newaxis]]
print A

```

```py
[[ 0.90483742  0.1       ]
 [ 0.81873075  0.2       ]
 [ 0.74081822  0.3       ]
 [ 0.67032005  0.4       ]
 [ 0.60653066  0.5       ]
 [ 0.54881164  0.6       ]
 [ 0.4965853   0.7       ]
 [ 0.44932896  0.8       ]
 [ 0.40656966  0.9       ]
 [ 0.36787944  1\.        ]]

```

求解最小二乘问题：

In [16]:

```py
c, resid, rank, sigma = linalg.lstsq(A, zi)

print c

```

```py
[ 4.87016856  2.19081311]

```

其中 `c` 的形状与 `zi` 一致，为最小二乘解，`resid` 为 `zi - A c` 每一列差值的二范数，`rank` 为矩阵 `A` 的秩，`sigma` 为矩阵 `A` 的奇异值。

查看拟合效果：

In [17]:

```py
xi2 = np.r_[0.1:1.0:100j]
yi2 = c[0]*np.exp(-xi2) + c[1]*xi2

plt.plot(xi,zi,'x',xi2,yi2)
plt.axis([0,1.1,3.0,5.5])
plt.xlabel('$x_i$')
plt.title('Data fitting with linalg.lstsq')
plt.show()

```

![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXAAAAEbCAYAAADDKt+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJxuEsO+7IAWlgIoLgiim1lrBpa36KNcuLrWt16rXLvZqe6+K94qttd5bba11a6XWn9T+1Cp11xoraFGWIAhUVtnDmrAkkO1z/zgnOKSTZJJMZnKS9/PxyCMzZ77nnM+ZTN7zne9ZxtwdERGJnox0FyAiIk2jABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgEuzmNlkM1tlZnvN7Atm9qKZXdbIZSwzsyktVWMC6/+RmT1cz+NXmNnbjVjeejM7K7z94/qW3cg6q83s6CbMl29mG5NRg7QuWekuQJrOzNYDfYFKoApYDvweeMgTOMDfzIYBa4Esd69uYhn/Bdzn7r8M7z8Xs/wrgKvc/YyYaY8BG939lppp7j62ietOCnf/Sc3tJD0nh597d7+zWcWlkJnNAEa4+9fTXYskRj3waHPgfHfvCgwFfgrcBDzayOVYM2oYSvDG0dY05zkRSQkFeBvh7vvcfQ4wHbjczMYAmNl5ZrbYzErMbIOZ3RYz29/C38Vmts/MTjWzEWb2VzPbaWY7zOwPZtYt3jrNbA1wNDAnHELJMbMCM7vKzI4FfgNMCpe9x8y+BXwF+Pdw2nPhcmKHHGaY2VNmNitc5jIzOylmnSeG27M3bPdHM/vvOur72MxODG9/NRyCGB3ev8rMno1Z5+NxnpO9ZjaRsEdtZneb2W4zW2tm5ybyd4ldtpkNC2u4LKxth5n9OKbtBDN7N3yutpjZL80su47l9jKzOeHf9T0zuyPRYR4zu8nMNoXbt9LMzgq350fA9PBvszhse4WZrQnbrjWzr4TTM83s5+E2rDGza8NtU6akkJ7sNsbd3wc2AaeHk/YDX3P3bsB5wDVm9oXwsZqhjW7u3sXd54f3ZwIDgNHAEGBGHesaAWwg/BTg7uUEYefuvhK4Gng3XHYPd38YeAK4K5xWU0ft4Z4LgCeBbsDzwK8AzCwHeBb4LdAjbPPFOPPXKADyw9tnAmvC3zX3C+LME/ucdHX3vxP0xk8FVgK9gJ+R+KeceLVNBkYBnwVuNbNjwumVwA3hOiaFj3+njuXeD+wD+gGXA5fVsa4jhOu6Fjg5/OR2DrDe3V8G7gRmh3+b8WaWB9wLnBu2nQQUhov6FsHr6QTgZOCSRNYvyaUAb5u2AD0B3P0td/8wvL0UmM0nIfZPwwTuvsbd33D3CnffCfxvTPvGqmsYoqHhibfd/eVwHP8PwPHh9IlAprv/0t2r3P1Z4L16lvMWn9R+OvCTmPtTwscTre1jd380rOn3wAAz69vAdtS1vNvd/ZC7fwAsIQhB3H2Ru7/n7tXu/jHwEHGeezPLBC4CbnP3g+6+AphVT+2xqoAOwBgzy3b3De6+NqbW2suoBsaZWa67F7l7zXDZl4H/dffN7r6HIPw17JRiCvC2aRCwGyAcFnnTzLabWTFBr7hXXTOaWT8zmx1+xC4BHq+vfQspirldCnQMP5oPBDbXaruRuoPjb8AZZtYfyAT+BEw2s6MIetiFdcwXz7aaG+5eGt7s3Ij54y6LYPvyAMxslJn9xcy2hs/9TOI/930IDkCIPbJkUyIrdvfVwHcJPlUVmdmTZjagjrYHCIbk/hXYEtZW82lhQK31b0hk/ZJcCvA2xsxOIQjwueGk/wf8GRjs7t0JxqVr/u7xPvLeSdBLGxsOu3ydpr9O4i2/OR+ztxJsW6yhdS0zDKtS4HrgLXffRxCe3wZix4u9jtup9gDBDuFPhc/9fxD/ud9BMNwyJGbakDjt4nL3J8Mjg44i2N67ah6K0/ZVdz8H6E8whFRzSORWgue+xtDa80rLU4BHnwGYWVczO59gXPjxmmETgl7iHncvN7MJBDsRa/5RdxB8RB4Rs7zOwAFgr5kNAn7YjNqKgMG1dsQVEez4bIp3gSozu87MssKx/FMamOct4Do+GS4pqHUfjuzBx3tOUqUzwbh2abgT+Jp4jdy9CngGmGFmuWHbr5PYGPiocKdlB+AQcJDgDRuCN7dhZlbzmuprwbH9eUAFweuipu1TwL+Z2SAz6wHcnMj6JbkU4NE3x8z2EnyE/RFwD3BlzOPfAf4rbHML8MeaB8KhgJnAvPDoignA7cCJQAkwB3iapv9jvgF8CGwzs+3htEeBT4dHWjwTZx6Psz4P6y0nGPu9CtgDfBX4C1BeTw1vEQTj3+q4f8Q64zwnp9ZXUwJqz1vffDcSvMHuJRj/nl3PvNcR7OTdRjD+/SQxz0N49M6lcebtQLAvYAdBL7o3wesGgiEmgF1mtoAgH75HMGy1i2AHb82bysPAKwRj+AsIXicaA08xS+QLHSw4YWQvwbtvhbtPqPV4PsEJHDU7Q5529zuSWqlIHGY2H/i1u89Kdy3pZGZ3AX3d/coGG7fM+ofR/BOgpJESPRPTgXx3311Pm7fc/cIk1CRSJwtOuf8I2EnQAx8LvJzWotIg3JnYAVhKMIz0DYJPJtKONOZU+oY+Hunjk6TCMQTjr3kEx3Vf4u5F9c/SJnUhGDYZSLBf4efu/nx6S9IYeKolOoSylmBMtAp4MDwhI/bxMwl2qmwiGC+7MeZ4URERaQGJ9sAnu/tWM+sDvGZmK9099jCsRcAQdy81s6kEh62NSnaxIiLyiYR64EfMEFxLY7+731NPm3XASbFj5mamj1ciIk3g7nGHqBs8jNDMOplZl/B2HsG1E5bWatMv5tjRCQRvDP+0w9Pd29XPbbfdlvYatM3aZm1ztLe3PokMofQDng3zOQt4wt1fNbOrw1B+kOBCNteYWSXBmW//ksByRUSkGRoMcHdfR3ixnVrTH4y5fT/B1dFERCRFdCZmC8rPz093CSmnbW4f2ts2t9btbfROzCavyMxTtS4RkbbCzPCm7sQUEZHWSQEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQkFuJmtN7MPzGyxmb1XR5v7zGyVmS0xs/HJLVNERGrLSrCdA/nuvjveg2Y2DfiUu480s1OBB4CJSapRRETiaMwQitXz2IXALAB3nw90N7N+zSlMRETql2iAO/C6mS0ws2/FeXwQsDHm/iZgcHOLExGRuiU6hDLZ3beaWR/gNTNb6e5v12pTu4fuzS9PRETqklCAu/vW8PcOM3sWmADEBvhmYEjM/cHhtCPMmDHj8O38/Hzy8/MbXbCISFtWUFBAQUFBQm3Nvf6Ospl1AjLdfZ+Z5QGvAre7+6sxbaYB17n7NDObCPzC3SfWWo43tC4RETmSmeHucfdBJtID7wc8a2Y17Z9w91fN7GoAd3/Q3V80s2lmtho4AFyZpNpFRKQODfbAk7Yi9cBFRBqtvh64zsQUEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUW06wKu9mlmFs6ioqkh3KSIiSdemA3zvob08sfQJTn3kVBZtXZTuckREkqpNB3j3jt155WuvcMOpNzD1ianc/PrNlFaUprssEZGkaNMBDmBmXH7C5Xzwrx+woWQD4x4Yx2trXkt3WSIizWbunpoVmXmq1lWfl1a9xDUvXMPpQ0/nfz7/P/TN65vukkRE6mRmuLvFe6zN98BrmzpyKh9+50MGdhnI2F+P5aGFD1Ht1ekuq9leeAGKi4+cVlwcTBeRtqndBThAXk4eP/vcz3j9std5rPAxJv92MoXbCtNdVrNMngz/8R+fhHhxcXB/8uT01iUiLSehADezTDNbbGZz4jyWb2Yl4eOLzew/k19m4yTaGz2u33HM/cZcrhp/FZ//w+f57svfpeRgSeoKTaLu3WHmzCC0168Pfs+cGUwXkbYp0R74DcByoK5B7LfcfXz4c0dySmu6xvRGMyyDb574TT78zofsL9/P6PtH8/iSx2kN4/WN1b07/PCHMHx48FvhLdK2NRjgZjYYmAY8AsQdSK9nelo0pTfau1NvHrnwEZ6d/iz3zr+XKY9NidywSnEx3H03rFsX/K79KURE2pYGj0Ixsz8BdwJdgRvd/YJaj58JPANsAjaHbZbHWU7Kj0JZvz7oja5bB8OGJT5fVXUVjy5+lFvfvJUvHfsl7jjrDnp16tVSZSZFzaeMmjeq2vdFJJqafBSKmZ0PbHf3xdTdy14EDHH344FfAn9uTrHJ0pzeaGZGJt8+6dusuHYF2ZnZjL5/NPfNv69Vn5I/b96RYV3zKWTevPTWJSItp94euJndCXwdqAQ6EvTCn3b3y+qZZx1wkrvvrjXdb7vttsP38/Pzyc/Pb1bxdUl2b/TD7R/yvVe+x6a9m7jnnHuYOnJq8osWEQEKCgooKCg4fP/222+vswee8Ik84VBJvCGUfgS9dDezCcBT7j4szvwpG0J54YVgh2VsWBcXB73R885r2jLdnTkfzeHGV2/k6B5Hc8859zCm75jkFCwiUof6hlAaG+A/cPcLzexqAHd/0MyuBa4h6KWXAt9397/Hmb9VnInZXOVV5Tzw/gPMfHsmF42+iNvzb6df537pLktE2qikBHgSimgTAV5jT9keZr49k8cKH+OGU2/g+5O+T15OXrrLEpE2RqfSt4AeuT34+Tk/571vvcfyncsZ9atRPLTwISqrK9Ndmoi0E+qBJ8mCLQu46fWb2Lx3M3ecdQcXj74Ys1Z1eLyIRJCGUFLE3Xlt7Wvc/PrNZGZkcudZd3L20WcryEWkyRTgKVbt1fzpwz9xy5u3MKjrIGaeNZPThpyW7rJEJII0Bp5iGZbB9LHTWX7tcr427mtc+vSlTHtiGgu3LEx3aS1Cl7IVSQ8FeAvKysjiqhOv4qPrPuK8kedx4ewL+eLsL0buGisN0aVsRdJDQygpVFZRxkMLH+KueXcxcfBEbj3zVk7of0K6y0qKmtD+4Q+DSxfoGiwiyaEx8FamrKKM3yz4DXe/czenDDqFW6bcwskDT053Wc3W1IuHiUjdNAbeyuRm5/K9Sd9jzb+t4ezhZ/OlP36JqU9MZe6Guekurcl0KVuR1FOAJ1ljdujlZudy/anXs/r61Vx07EVc/ufLOfOxM3l59cuR+kKJ2IuFDRv2ybXYFeIiLUtDKEnWnCshVlZXMnvZbO6adxdZGVncNPkmLvn0JWRlZKWm+CZqiYuHiUhAY+Ap1twdeu7Oi6te5KfzfsrmvZv5waQfcOX4K+mU3anlihaRVkkBngbJ2qH37sZ3ufudu5m7YS5Xn3Q11024Tlc/FGlHtBMzxZK5Q2/SkEk8M/0Z5n5jLjtLdzL6/tFc9dxVLC1amryCRSSS1ANPspb+bsqdpTt5cMGD3P/+/YzpO4YbTr2BaSOnkWF6LxZpizSEkkKp2qFXXlXOH5f9kXvn30vxwWKum3AdV5xwBd076uwZkbZEAd6GuTt/3/R37nvvPl5e/TLTx0zn2lOuZVy/cekuTUSSQAHeTmzdt5WHFz3Mgwsf5OgeR3PNyddw8eiL6ZDVId2liUgTKcDbmYqqCuZ8NIcHFjzAkm1LuPz4y/n2Sd9mZK+R6S5NRBpJAd6Ord69mocXPsxjSx5jTJ8xfPPEb3LR6IvomNUx3aWJSAIU4EJ5VTnPrXyORxY/wsItC7l07KV8Y/w3GD9gfLpLE5F6KMDlCOuL1zOrcBa/K/wdPXJ7cMXxV/CVcV+hT16fdJcmIrUowCWuaq/mzXVvMmvJLJ7/x/OcOexMLjvuMs4fdb52fIq0EgpwadC+Q/t4esXTPP7B4xRuK+Ti0Rfz1XFf5YyjztBJQiJppACXRtlYspEnlz3JE0ufYE/ZHqaPmc6l4y5lfP/xmMV9HYlIC1GAS5Mt276M2ctm8+SyJ8m0TKaPmc6Xx3yZsX3HKsxFUkABLs3m7ry/5X2e+vApnvrwKfJy8rhk9CVc8ulLOK7fcQrzdkDXfU8PBbgkVbVX897m93h6+dM8veJpMjMy+dKxX+Ki0RcxYdAEjZm3US19oTaJTwEuLcbdKdxWyDMrnuGZlc+wp2wPF4y6gC8c+wXOGn6WThhqY5r7ZSXSeApwSZlVu1bx3D+e488r/8zS7Uv57PDPcsGoC5g2cpq+iKKNSNaXlUhiFOCSFjtLd/LiqheZ89EcXl/7OiN7jmTayGlM/dRUTh54MpkZmekuURpJPfDUU4BL2pVXlTN3w1xeWvUSL61+iaIDRZx99Nl8fsTnOWfEOQzsMjDdJUoDNAaeHgpwaXU2lGzg1TWv8sqaV3hj7RsM6jqIs4efzedGfI4pR02hc07ndJcotaTrKJT2fvSLAlxatarqKhZuXchra17j9XWv8/7m9zm+//GcNewsPjP8M0waPInc7Nx0lylp0t57/gpwiZTSilLe2fgOf133VwrWF/BB0QecOOBEphw1hSlHTWHS4El06dAl3WVKCrXnsXcFuETa/vL9vLvxXd76+C3+9vHfWLR1Ecf0PoYzhp7BaUNOY/KQyQzqOijdZUoLa69HvzQ7wM0sE1gAbHL3C+I8fh8wFSgFrnD3xXHaKMAlKQ5VHmLh1oW8/fHbvLPpHd7Z+A65WblMGjKJiYMmMnHwRE7of4KGXdoQ9cCbF+DfB04Curj7hbUemwZc5+7TzOxU4F53nxhnGQpwaRHuzqrdq5i/aT7vbnqX+Zvns2LHCkb3Gc0pA0/h5IEnc/LAkxnTZwzZmdnpLlcaSWPgzQhwMxsMPAbMBL5fuwduZr8B3nT3P4b3VwJnuntRrXYKcEmZsooyCrcVsmDLAt7f8j4LtixgffF6xvQdw4n9T2T8gPGc0P8ExvUdR15OXrrLlXroKJTmBfifgDuBrsCNcQJ8DvATd38nvP86cJO7L6zVTgEuabW/fD9Lti1h0dZFFG4rZPG2xazcuZIh3YZwXL/jOK7vcYztO5Zx/cYxvPtwnWgkrUJ9AZ7VwIznA9vdfbGZ5dfXtNZ9JbW0Op1zOjN56GQmD518eFpFVQX/2PUPlmxbwrLty/ht4W9ZWrSU7Qe2c0zvYxjTZwyje49mdJ/RHNv7WEb0GKFvK5JWo94AB04DLgzHuTsCXc3s9+5+WUybzcCQmPuDw2n/ZMaMGYdv5+fnk5+f34SSRZInOzObsX3HMrbv2COm7y/fz4odK1i+Yzkrdq7gscLHWLlzJRtKNjCo6yBG9RrFyJ4jGdVrFCN6jGBEzxEM6z6MnMycNG2JtCR3T9klkwsKCigoKEiobcKHEZrZmcQfQondiTkR+IV2YkpbVVFVwdo9a1m1exWrdq3io10fsWbPGtbsWcOmvZvo37k/R/c4muHdh3NUt6MY1n0YR3U/iqHdhjK462AFfCtUWV1J0f4ituzbwuZ9m9m0d9Phnw0lG9i4dyPZGdl8dP1HaakvKceBhwH+A3e/0MyuBnD3B8PHfgWcCxwArnT3RXHmV4BLm1ZRVcHGvRtZu2ct6/as4+OSj1lfvJ4NJRvYULKBLfu20DO3J0O6DWFQl0HBT9dBDOg8gAFdBjCg8wD6de5Hn059NP7eTNVeTfHBYrYf2E7R/qLg94EiivYXsW3/Nrbu3xr87NvKjtId9MrtxaCugxjYZSCDuww+/Dca2m0oQ7oNYXDXwWm7NLJO5BFpBaqqqyg6UHS4d7d572Y279vMln1bDofKtv3b2FO2h565Pemb15c+eX3o06kPvTv1pnen3vTK7UWvTr3omduTHh170L1jd7p37E63jt3IzcptU9+M5O4cqjrE3kN7KTlYQvHB4sM/ew7uYU/ZHnaX7WZ32W52le1iV9kudpbuZGfpTnaX7aZzTmf65vWlb15f+uX1C34692NA5wH079yf/p37M7DLQPrm9W3Vh5cqwEUipLK6kh0HdrCjdMfh3ztLd7KrdBc7SnccEV4lh0rYU7aHkkMlVFZX0rVDV7rkdKFLhy50yelC55zOdM7pTF5OHp2yOpGbnUtuVi4dszqSm51Lh8wO5GTmkJOZQ3ZmNtkZ2WRnZpOVkUVWRhaZlkmGZZBhGYffHAzDcdwdx6mqrqLaq6nyKiqrK6msrqSiqoKK6grKq8opryrnUOUhDlUd4mDlQcoqyiirLDv8+0DFAQ6UH2B/+f7DP/vK97Hv0D4AunXsRrcO3ejWsRvdO3Y//MbVM7fn4Tey3p1606tTL3rl9qJPXh965fZq1aHcGApwkXagvKqckoMlh8NvX/m+w8F4oOIAZRVBWB6sPHg4SMuryjlUdYjyqnIqqisOB29VdRUV1RVUe3UQztVVAIeD28wwDDMj0zLJzAiCPjsjm8yMzMNvBNkZ2XTI7ECHrA50yOxAbnbw5tExqyOdsjuRm5VLXk4eedl55OXkHfGm07VDVx3xgwJcRKTRWssJRPUFuL59VkQkjsmTg1P2i4uD+zWn8E+eXP98qaQeuIhIHVrDRbQ0hCIi0kTpvoythlCkTXrhhU8+3tYoLg6miyRDcXHQ8163Lvhd+/WWbgpwiawojFFKdMVetnbYsOB37OutNdAQikRaaxijlLYpCkehKMAl8tI9RinSkjQGLm1Wax+jFGlJCnCJrHSOUaZrB6p23EosBbhE1rx5R455d+8e3J83r+XXna4dqNpxK7E0Bi7SROnagaodt+2LdmKKtJB07UDVjtv2QzsxRVpAunagaset1FCAizRBunagRuHkEkkdDaGINEG6TvJoLSeXSOpoDFxEJKI0Bi4i0gYpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRDUY4GbW0czmm1mhmS03s5/EaZNvZiVmtjj8+c+WKVdERGpkNdTA3Q+a2WfcvdTMsoC5Zna6u8+t1fQtd7+wZcoUEZHaEhpCcffS8GYOkAnsjtMs7ne2iYhIy0gowM0sw8wKgSLgTXdfXquJA6eZ2RIze9HMPp3sQkVE5EiJ9sCr3f0EYDAwxczyazVZBAxx9+OBXwJ/TmqVIiLyTxocA4/l7iVm9gJwMlAQM31fzO2XzOzXZtbT3Y8YapkxY8bh2/n5+eTn5zetahGRNqqgoICCgoKE2pq719/ArDdQ6e7FZpYLvALc7u5vxLTpB2x3dzezCcBT7j6s1nK8oXWJiMiRzAx3j7uPMZEe+ABglpllEAy5PO7ub5jZ1QDu/iBwCXCNmVUCpcC/JKd0ERGpS4M98KStSD1wEZFGq68HrjMxRUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIqjfAzayjmc03s0IzW25mP6mj3X1mtsrMlpjZ+JYpVUREYtUb4O5+EPiMu58AHAd8xsxOj21jZtOAT7n7SODbwAMtVWzUFBQUpLuElNM2tw/tbZtb6/Y2OITi7qXhzRwgE9hdq8mFwKyw7Xygu5n1S2aRUdVa/+gtSdvcPrS3bW6t29tggJtZhpkVAkXAm+6+vFaTQcDGmPubgMHJK1FEROJJpAdeHQ6hDAammFl+nGZWe7Yk1CYiIvUw98Sz1sxuAcrc/ecx034DFLj77PD+SuBMdy+qNa9CXUSkCdy9dicZgKz6ZjKz3kCluxebWS7wOeD2Ws2eB64DZpvZRKC4dnjXV4CIiDRNvQEODABmmVkGwXDL4+7+hpldDeDuD7r7i2Y2zcxWAweAK1u2ZBERgUYOoYiISOuR9DMxzexcM1sZnthzUx1t2tSJPw1ts5l9NdzWD8xsnpkdl446kyWRv3HY7hQzqzSzi1JZX0tI8HWdb2aLzWyZmRWkuMSkS+B13dvMXg5P9FtmZlekocykMbPfmlmRmS2tp03ryi53T9oPwXHiq4FhQDZQCIyu1WYa8GJ4+1Tg78msIdU/CW7zJKBbePvcKG9zItsb0+6vwF+Ai9Nddwr+xt2BD4HB4f3e6a47Bds8A/hJzfYCu4CsdNfejG0+AxgPLK3j8VaXXcnugU8AVrv7enevAGYDX6jVpq2d+NPgNrv7u+5eEt6dT7SPk0/kbwxwPfD/gR2pLK6FJLLNXwGedvdNAO6+M8U1Jlsi27wV6Bre7grscvfKFNaYVO7+NrCnniatLruSHeDxTuoZlECbKAdaItsc6yrgxRatqGU1uL1mNojgn73msgpR39GSyN94JNDTzN40swVm9vWUVdcyEtnmh4ExZrYFWALckKLa0qXVZVdDR6E0VqL/qG3pxJ+EazezzwDfACa3XDktLpHt/QVws7u7mRn//PeOmkS2ORs4Efgs0Al418z+7u6rWrSylpPINv8YKHT3fDMbAbxmZse7+74Wri2dWlV2JTvANwNDYu4PIXiXqq/N4HBaVCWyzYQ7Lh8GznX3+j6mtXaJbO9JBOcFQDA2OtXMKtz9+dSUmHSJbPNGYKe7lwFlZvY34HggqgGeyDafBswEcPc1ZrYOOAZYkJIKU6/VZVeyh1AWACPNbJiZ5QDTCU70ifU+xBxCAAABwklEQVQ8cBlAfSf+REiD22xmQ4FngK+5++o01JhMDW6vux/t7sPdfTjBOPg1EQ5vSOx1/Rxwupllmlkngp1cta8bFCWJbPNK4GyAcCz4GGBtSqtMrVaXXUntgbt7pZldB7xCsBf7UXdf0ZZP/Elkm4FbgR7AA2GvtMLdJ6Sr5uZIcHvblARf1yvN7GXgA6AaeNj/+cJvkZHg3/lO4HdmtoSgM/jv7l77aqWRYWZPAmcCvc1sI3AbwdBYq80uncgjIhJR+ko1EZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElHJvhaKSCSYWSbB6eFHE1zHZAJwj7u35VPBpY1RD1zaq+OBpwmu3ZEB/Ing+tYikaEAl3bJ3Re5+yGCb0sqcPcCdy8zs8+nuzaRRCnApV0Kv6+zNzDW3deZ2ekA7v5KmksTSZguZiXtkpndAhQBQ4GFwHbgEHC6u/8inbWJJEo7MaVdcvf/rj0t/FaZ4jSUI9IkGkIR+cSJKMAlQjSEIiISUeqBi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRNT/AXU6OlSJKzNiAAAAAElFTkSuQmCC)

### 广义逆

`linalg.pinv` 或 `linalg.pinv2` 可以用来求广义逆，其区别在于前者使用求最小二乘解的算法，后者使用求奇异值的算法求解。

## 矩阵分解

### 特征值和特征向量

#### 问题描述

对于给定的 $N \times N$ 矩阵 $\mathbf A$，特征值和特征向量问题相当与寻找标量 $\lambda$ 和对应的向量 $\mathbf v$ 使得： $$ \mathbf{Av} = \lambda \mathbf{v} $$

矩阵的 $N$ 个特征值（可能相同）可以通过计算特征方程的根得到： $$ \left|\mathbf{A} - \lambda \mathbf{I}\right| = 0 $$

然后利用这些特征值求（归一化的）特征向量。

#### 问题求解

*   `linalg.eig(A)`
    *   返回矩阵的特征值与特征向量
*   `linalg.eigvals(A)`
    *   返回矩阵的特征值
*   `linalg.eig(A, B)`
    *   求解 $\mathbf{Av} = \lambda\mathbf{Bv}$ 的问题

#### 例子

矩阵为 $$ \begin{split}\mathbf{A}=\left[\begin{array}{ccc} 1 & 5 & 2\\ 2 & 4 & 1\\ 3 & 6 & 2\end{array}\right].\end{split} $$

特征多项式为： $$ \begin{eqnarray*} \left|\mathbf{A}-\lambda\mathbf{I}\right| & = & \left(1-\lambda\right)\left[\left(4-\lambda\right)\left(2-\lambda\right)-6\right]-\\ & & 5\left[2\left(2-\lambda\right)-3\right]+2\left[12-3\left(4-\lambda\right)\right]\\ & = & -\lambda^{3}+7\lambda^{2}+8\lambda-3.\end{eqnarray*} $$

特征根为： $$ \begin{eqnarray*} \lambda_{1} & = & 7.9579\\ \lambda_{2} & = & -1.2577\\ \lambda_{3} & = & 0.2997.\end{eqnarray*} $$

In [18]:

```py
A = np.array([[1, 5, 2], 
              [2, 4, 1],
              [3, 6, 2]])

la, v = linalg.eig(A)

print la

# 验证是否归一化
print np.sum(abs(v**2),axis=0)

# 第一个特征值
l1 = la[0]
# 对应的特征向量
v1 = v[:, 0].T

# 验证是否为特征值和特征向量对
print linalg.norm(A.dot(v1)-l1*v1)

```

```py
[ 7.95791620+0.j -1.25766471+0.j  0.29974850+0.j]
[ 1\.  1\.  1.]
3.23301824835e-15

```

### 奇异值分解

#### 问题描述

$M \times N$ 矩阵 $\mathbf A$ 的奇异值分解为： $$ \mathbf{A=U}\boldsymbol{\Sigma}\mathbf{V}^{H} $$

其中 $\boldsymbol{\Sigma}, (M \times N)$ 只有对角线上的元素不为 0，$\mathbf U, (M \times M)$ 和 $\mathbf V, (N \times N)$ 为正交矩阵。

其具体原理可以查看维基百科： [https://en.wikipedia.org/wiki/Singular_value_decomposition](https://en.wikipedia.org/wiki/Singular_value_decomposition)

#### 问题求解

*   `U,s,Vh = linalg.svd(A)`
    *   返回 $U$ 矩阵，奇异值 $s$，$V^H$ 矩阵
*   `Sig = linalg.diagsvd(s,M,N)`
    *   从奇异值恢复 $\boldsymbol{\Sigma}$ 矩阵

#### 例子

奇异值分解：

In [19]:

```py
A = np.array([[1,2,3],[4,5,6]])

U, s, Vh = linalg.svd(A)

```

$\boldsymbol{\Sigma}$ 矩阵：

In [20]:

```py
M, N = A.shape
Sig = linalg.diagsvd(s,M,N)

print Sig

```

```py
[[ 9.508032    0\.          0\.        ]
 [ 0\.          0.77286964  0\.        ]]

```

检查正确性：

In [21]:

```py
print A
print U.dot(Sig.dot(Vh))

```

```py
[[1 2 3]
 [4 5 6]]
[[ 1\.  2\.  3.]
 [ 4\.  5\.  6.]]

```

### LU 分解

$M \times N$ 矩阵 $\mathbf A$ 的 `LU` 分解为： $$ \mathbf{A}=\mathbf{P}\,\mathbf{L}\,\mathbf{U} $$

$\mathbf P$ 是 $M \times M$ 的单位矩阵的一个排列，$\mathbf L$ 是下三角阵，$\mathbf U$ 是上三角阵。

可以使用 `linalg.lu` 进行 LU 分解的求解：

具体原理可以查看维基百科： [https://en.wikipedia.org/wiki/LU_decomposition](https://en.wikipedia.org/wiki/LU_decomposition)

In [22]:

```py
A = np.array([[1,2,3],[4,5,6]])

P, L, U = linalg.lu(A)

print P
print L
print U

print P.dot(L).dot(U)

```

```py
[[ 0\.  1.]
 [ 1\.  0.]]
[[ 1\.    0\.  ]
 [ 0.25  1\.  ]]
[[ 4\.    5\.    6\.  ]
 [ 0\.    0.75  1.5 ]]
[[ 1\.  2\.  3.]
 [ 4\.  5\.  6.]]

```

### Cholesky 分解

`Cholesky` 分解是一种特殊的 `LU` 分解，此时要求 $\mathbf A$ 为 Hermitian 正定矩阵 （$\mathbf A = \mathbf{A^H}$）。

此时有： $$ \begin{eqnarray*} \mathbf{A} & = & \mathbf{U}^{H}\mathbf{U}\\ \mathbf{A} & = & \mathbf{L}\mathbf{L}^{H}\end{eqnarray*} $$ 即 $$ \mathbf{L}=\mathbf{U}^{H}. $$

可以用 `linalg.cholesky` 求解。

### QR 分解

$M×N$ 矩阵 $\mathbf A$ 的 `QR` 分解为： $$ \mathbf{A=QR} $$

$\mathbf R$ 为上三角形矩阵，$\mathbf Q$ 是正交矩阵。

维基链接： [https://en.wikipedia.org/wiki/QR_decomposition](https://en.wikipedia.org/wiki/QR_decomposition)

可以用 `linalg.qr` 求解。

### Schur 分解

对于 $N\times N$ 方阵 $\mathbf A$, `Schur` 分解要求找到满足下式的矩阵： $$ \mathbf{A=ZTZ^H} $$

其中 $\mathbf Z$ 是正交矩阵，$\mathbf T$ 是一个上三角矩阵。

维基链接： [https://en.wikipedia.org/wiki/Schur_decomposition](https://en.wikipedia.org/wiki/Schur_decomposition)

In [23]:

```py
A = np.mat('[1 3 2; 1 4 5; 2 3 6]')

print A

T, Z = linalg.schur(A)

print T, Z

print Z.dot(T).dot(Z.T)

```

```py
[[1 3 2]
 [1 4 5]
 [2 3 6]]
[[ 9.90012467  1.78947961 -0.65498528]
 [ 0\.          0.54993766 -1.57754789]
 [ 0\.          0.51260928  0.54993766]] [[ 0.36702395 -0.85002495 -0.37782404]
 [ 0.63681656 -0.06646488  0.76814522]
 [ 0.67805463  0.52253231 -0.51691576]]
[[ 1\.  3\.  2.]
 [ 1\.  4\.  5.]
 [ 2\.  3\.  6.]]

```

## 矩阵函数

考虑函数 $f(x)$ 的泰勒展开： $$ f\left(x\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}x^{k} $$

对于方阵，矩阵函数可以定义如下： $$ f\left(\mathbf{A}\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}\mathbf{A}^{k} $$

这也是计算矩阵函数的最好的方式。

### 指数和对数函数

#### 指数

指数可以定义如下： $$ e^{\mathbf{A}}=\sum_{k=0}^{\infty}\frac{1}{k!}\mathbf{A}^{k} $$

`linalg.expm3` 使用的是泰勒展开的方法计算结果：

In [24]:

```py
A = np.array([[1, 2], [3, 4]])

print linalg.expm3(A)

```

```py
[[  51.96890355   74.73648784]
 [ 112.10473176  164.07363531]]

```

另一种方法先计算 A 的特征值分解： $$ \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} $$

然后有（正交矩阵和对角阵的性质）： $$ e^{\mathbf{A}}=\mathbf{V}e^{\boldsymbol{\Lambda}}\mathbf{V}^{-1} $$

`linalg.expm2` 使用的就是这种方法：

In [25]:

```py
print linalg.expm2(A)

```

```py
[[  51.9689562    74.73656457]
 [ 112.10484685  164.07380305]]

```

最优的方法是用 [`Padé` 近似](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant) 实现，`Padé` 近似往往比截断的泰勒级数准确，而且当泰勒级数不收敛时，`Padé` 近似往往仍可行，所以多用于在计算机数学中。

`linalg.expm` 使用的就是这种方法：

In [26]:

```py
print linalg.expm(A)

```

```py
[[  51.9689562    74.73656457]
 [ 112.10484685  164.07380305]]

```

#### 对数

指数的逆运算，可以用 `linalg.logm` 实现：

In [27]:

```py
print A
print linalg.logm(linalg.expm(A))

```

```py
[[1 2]
 [3 4]]
[[ 1\.  2.]
 [ 3\.  4.]]

```

### 三角函数

根据欧拉公式，其定义为： $$ \begin{eqnarray*} \sin\left(\mathbf{A}\right) & = & \frac{e^{j\mathbf{A}}-e^{-j\mathbf{A}}}{2j}\\ \cos\left(\mathbf{A}\right) & = & \frac{e^{j\mathbf{A}}+e^{-j\mathbf{A}}}{2}.\end{eqnarray*} $$

正切函数定义为： $$ \tan\left(x\right)=\frac{\sin\left(x\right)}{\cos\left(x\right)}=\left[\cos\left(x\right)\right]^{-1}\sin\left(x\right) $$

因此矩阵的正切函数定义为： $$ \left[\cos\left(\mathbf{A}\right)\right]^{-1}\sin\left(\mathbf{A}\right). $$

具体实现：

*   `linalg.sinm`
*   `linalg.cosm`
*   `linalg.tanm`

### 双曲三角函数

\begin{eqnarray*} \sinh\left(\mathbf{A}\right) & = & \frac{e^{\mathbf{A}}-e^{-\mathbf{A}}}{2}\\ \cosh\left(\mathbf{A}\right) & = & \frac{e^{\mathbf{A}}+e^{-\mathbf{A}}}{2}\\ \tanh\left(\mathbf{A}\right) & = & \left[\cosh\left(\mathbf{A}\right)\right]^{-1}\sinh\left(\mathbf{A}\right).\end{eqnarray*}

具体实现：

*   `linalg.sinhm`
*   `linalg.coshm`
*   `linalg.tanhm`

## 特殊矩阵

`Scipy` 提供了一些特殊矩阵的实现，具体可以参考：

[http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#special-matrices](http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#special-matrices)