| jupytext |
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| kernelspec |
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| translation |
|
(linear_algebra)=
<div id="qe-notebook-header" align="right" style="text-align:right;">
<a href="https://quantecon.org/" title="quantecon.org">
<img style="width:250px;display:inline;" width="250px" src="https://assets.quantecon.org/img/qe-menubar-logo.svg" alt="QuantEcon">
</a>
</div>
:depth: 2
线性代数是经济学家最值得投入学习的应用数学分支之一。
例如,经济学和金融学中的许多应用问题都需要求解线性方程组,比如
或者更一般地,
:label: la_se
\begin{aligned}
y_1 = a_{11} x_1 + a_{12} x_2 + \cdots + a_{1k} x_k \\
\vdots \\
y_n = a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nk} x_k
\end{aligned}
这里的目标是在已知
在研究这类问题时,我们需要考虑以下几个基本问题:
- 这个方程组是否有解?
- 如果有解,解是唯一的吗?如果有多个解,这意味着什么?
- 如果没有精确解,是否可以找到一个最优的近似解?
- 如果解存在,有什么有效的方法可以求出它?
这些都是线性代数所要解决的问题。
本讲将介绍线性代数和矩阵运算的基础知识。
我们将同时关注理论基础和实际计算。
虽然本讲与我们在之前的讲座中已经介绍过NumPy数组的基本操作有所重合,但这里我们将从更理论的角度来探讨线性代数。
请注意,本讲比大多数讲座更偏重理论,其中包含的背景知识将在我们后续的应用中逐步用到。
让我们从一些导入开始:
import matplotlib.pyplot as plt
import matplotlib as mpl
FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf"
mpl.font_manager.fontManager.addfont(FONTPATH)
plt.rcParams['font.family'] = ['Source Han Serif SC']
plt.rcParams["figure.figsize"] = (11, 5) #设置默认图形大小
import numpy as np
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from scipy.linalg import inv, solve, det, eig
长度为$n$的向量就是一个由$n$个数字组成的序列(或数组、元组),我们将其写作$x = (x_1, \ldots, x_n)$或$x = [x_1, \ldots, x_n]$。
我们可以根据需要将这些序列横向或纵向书写。
(在后面讨论矩阵运算时,区分横向和纵向表示将变得很重要)
我们用$\mathbb R^n$表示所有$n$维向量的集合。
以$\mathbb R^2$为例,它代表二维平面,其中的每个向量都对应平面上的一个点。
在传统上,向量在视觉上被表示为从原点指向该点的箭头。
下图以这种方式展示了三个向量
fig, ax = plt.subplots(figsize=(10, 8))
# 设置通过原点的轴
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
ax.set(xlim=(-5, 5), ylim=(-5, 5))
ax.grid()
vecs = ((2, 4), (-3, 3), (-4, -3.5))
for v in vecs:
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(facecolor='blue',
shrink=0,
alpha=0.7,
width=0.5))
ax.text(1.1 * v[0], 1.1 * v[1], str(v))
plt.show()
向量最常见的两种运算是加法和标量乘法,下面我们来介绍。
根据定义,当我们对两个向量进行加法运算时,我们是按元素逐个相加
标量乘法是一种将数字
下图展示了标量乘法
---
tags: [output_scroll]
---
fig, ax = plt.subplots(figsize=(10, 8))
# Set the axes through the origin
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
ax.set(xlim=(-5, 5), ylim=(-5, 5))
x = (2, 2)
ax.annotate('', xy=x, xytext=(0, 0),
arrowprops=dict(facecolor='blue',
shrink=0,
alpha=1,
width=0.5))
ax.text(x[0] + 0.4, x[1] - 0.2, '$x$', fontsize='16')
scalars = (-2, 2)
x = np.array(x)
for s in scalars:
v = s * x
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(facecolor='red',
shrink=0,
alpha=0.5,
width=0.5))
ax.text(v[0] + 0.4, v[1] - 0.2, f'${s} x$', fontsize='16')
plt.show()
在Python中,向量可以用列表或元组来表示,例如 x = (2, 4, 6),但更常见的是用NumPy数组来表示。
NumPy数组的一个优点是标量乘法和加法具有非常自然的语法
x = np.ones(3) # 创建一个包含三个1的向量
y = np.array((2, 4, 6)) # 将元组(2, 4, 6)转换为数组
x + y
4 * x
向量
如果两个向量的内积为零,则称这两个向量正交。
向量
表达式
继续前面的例子,内积和范数可以按如下方式计算
np.sum(x * y) # x和y的内积,方法1
x @ y # x和y的内积,方法2(推荐)
推荐使用@运算符,因为它使用了经过优化的BLAS库来实现乘加融合运算,相比分别进行乘法和求和运算,具有更好的性能和数值精度。
np.sqrt(np.sum(x**2)) # x的范数,方法1
np.sqrt(x @ x) # x的范数,方法2(推荐)
np.linalg.norm(x) # x的范数,方法3
给定$\mathbb R ^n$中的一组向量$A := {a_1, \ldots, a_k}$,我们自然会思考通过执行线性运算可以创建哪些新向量。
以这种方式创建的新向量被称为$A$的线性组合。
具体来说,如果$y \in \mathbb R ^n$是$A := {a_1, \ldots, a_k}$的线性组合,那么
在这种情况下,值$\beta_1, \ldots, \beta_k$被称为线性组合的系数。
下图显示了$\mathbb R ^3$中$A = {a_1, a_2}$的张成空间。
张成空间是一个通过这两个点和原点的二维平面。
(la_3dvec)=
---
tags: [output_scroll]
---
ax = plt.figure(figsize=(10, 8)).add_subplot(projection='3d')
x_min, x_max = -5, 5
y_min, y_max = -5, 5
α, β = 0.2, 0.1
ax.set(xlim=(x_min, x_max), ylim=(x_min, x_max), zlim=(x_min, x_max),
xticks=(0,), yticks=(0,), zticks=(0,))
gs = 3
z = np.linspace(x_min, x_max, gs)
x = np.zeros(gs)
y = np.zeros(gs)
ax.plot(x, y, z, 'k-', lw=2, alpha=0.5)
ax.plot(z, x, y, 'k-', lw=2, alpha=0.5)
ax.plot(y, z, x, 'k-', lw=2, alpha=0.5)
# 固定线性函数,生成平面
def f(x, y):
return α * x + β * y
# 向量位置,按坐标
x_coords = np.array((3, 3))
y_coords = np.array((4, -4))
z = f(x_coords, y_coords)
for i in (0, 1):
ax.text(x_coords[i], y_coords[i], z[i], f'$a_{i+1}$', fontsize=14)
# 连接向量的线
for i in (0, 1):
x = (0, x_coords[i])
y = (0, y_coords[i])
z = (0, f(x_coords[i], y_coords[i]))
ax.plot(x, y, z, 'b-', lw=1.5, alpha=0.6)
# 绘制平面
grid_size = 20
xr2 = np.linspace(x_min, x_max, grid_size)
yr2 = np.linspace(y_min, y_max, grid_size)
x2, y2 = np.meshgrid(xr2, yr2)
z2 = f(x2, y2)
ax.plot_surface(x2, y2, z2, rstride=1, cstride=1, cmap=cm.jet,
linewidth=0, antialiased=True, alpha=0.2)
plt.show()
如果
如果
那么
现在考虑
如果
因此
(la_li)=
在实际应用中,我们经常需要找到一组能够张成较大空间的向量。这样做的好处是,我们可以用这少数几个向量的线性组合来表示更多的向量。
要使一组向量能张成较大的空间,这组向量需要满足一个重要条件 -- 线性无关。
具体来说,在$\mathbb R ^n$中的一组向量$A := {a_1, \ldots, a_k}$被称为
- 线性相关的,如果$A$的某个真子集与$A$具有相同的张成空间。
- 线性无关的,如果它不是线性相关的。
换句话说,一组向量如果没有一个向量对张成空间来说是多余的,则该组向量是线性无关的,否则就是线性相关的。
为了说明这个概念,让我们回到{ref}图示<la_3dvec>中的例子,其中$\mathbb R ^3$中的两个向量${a_1, a_2}$张成了一个经过原点的平面。
如果我们添加第三个向量$a_3$,构成集合${a_1, a_2, a_3}$,那么这个集合将是
- 线性相关的,如果$a_3$落在该平面上
- 线性无关的,否则
作为该概念的另一个例证,由于$\mathbb R ^n$可以由$n$个向量张成(参见上面关于标准基向量的讨论),因此在$\mathbb R ^n$中任何超过$n$个向量(即$m > n$个向量)的集合必定是线性相关的。
以下陈述与$A := {a_1, \ldots, a_k} \subset \mathbb R ^n$的线性无关性是等价的
- 集合中的任何向量都不能写成其他向量的线性组合。
- 如果对于标量$\beta_1, \ldots, \beta_k$,有$\beta_1 a_1 + \cdots + \beta_k a_k = 0$,那么$\beta_1 = \cdots = \beta_k = 0$。
(第一个表达式中的零指的是$\mathbb R ^n$的原点)
(la_unique_reps)=
线性独立向量集合的一个重要性质是唯一表示性:其张成空间中的每个向量都可以用唯一的一组系数来表示。
具体来说,假设$A := {a_1, \ldots, a_k} \subset \mathbb R ^n$是线性独立的,且
那么不存在其他系数序列$\gamma_1, \ldots, \gamma_k$能够得到相同的向量$y$。
事实上,如果我们同时有$y = \gamma_1 a_1 + \cdots \gamma_k a_k$,那么
由线性无关性可知,对所有的$i$都有$\gamma_i = \beta_i$。
矩阵是组织数据用于线性运算的一种简洁方式。
一个
通常,矩阵中的数字表示线性方程组中的系数,正如本讲开始时所讨论的。
显然,当
在前一种情况下,$A$ 被称为行向量,而在后一种情况下被称为列向量。
如果
将每个
如果
对于方阵
如果只有主对角线上的元素非零,则
如果除了是对角矩阵外,主对角线上的每个元素都等于1,则
与向量一样,矩阵也定义了许多代数运算。
标量乘法和加法是向量情况的直接推广:
以及
在后一种情况下,矩阵必须具有相同的形状才能使定义有意义。
我们还有一个相乘两个矩阵的约定。
矩阵乘法的规则推广了上面讨论的内积的概念。
并且其设计使得乘法能够很好地配合基本线性运算。
如果$A$和$B$是两个矩阵,那么它们的乘积$A B$的第$i,j$个元素是由$A$的第$i$行与$B$的第$j$列的内积得到的。
有许多教程可以帮助你理解这个运算,比如这个,或者维基百科页面上的讨论。
如果$A$是$n \times k$矩阵,$B$是$j \times m$矩阵,那么要使$A$和$B$可以相乘,我们需要$k = j$,且得到的矩阵$A B$是$n \times m$的。
作为可能最重要的特例,考虑将$n \times k$矩阵$A$与$k \times 1$列向量$x$相乘。
根据前面的规则,这会得到一个$n \times 1$列向量
:label: la_atx
A x =
\begin{bmatrix}
a_{11} & \cdots & a_{1k} \\
\vdots & \vdots & \vdots \\
a_{n1} & \cdots & a_{nk}
\end{bmatrix}
\begin{bmatrix}
x_{1} \\
\vdots \\
x_{k}
\end{bmatrix} :=
\begin{bmatrix}
a_{11} x_1 + \cdots + a_{1k} x_k \\
\vdots \\
a_{n1} x_1 + \cdots + a_{nk} x_k
\end{bmatrix}
矩阵乘法不满足交换律,也就是说,$A B$ 和 $B A$ 的结果通常是不同的。
另一个重要的特殊情况是单位矩阵。
你应该验证如果
如果
NumPy数组也被用作矩阵,并且对所有标准矩阵运算都有快速、高效的函数和方法1。
你可以通过元组的元组(或列表的列表)手动创建它们,如下所示
A = ((1, 2),
(3, 4))
type(A)
A = np.array(A)
type(A)
A.shape
shape 属性是一个表示行数和列数的元组 --- 更多讨论请参见这里。
要获取 A 的转置,可以使用 A.transpose() 或更简单的 A.T。
有许多便捷的函数可用于创建常见矩阵(零矩阵、全1矩阵等)--- 参见这里。
由于运算默认按元素执行,标量乘法和加法的语法具有非常自然
A = np.identity(3)
B = np.ones((3, 3))
2 * A
A + B
矩阵相乘我们使用@符号。
具体来说,A @ B表示矩阵乘法,而A * B表示元素逐个相乘。
更多讨论请参见这里。
(la_linear_map)=
每个$n \times k$矩阵$A$都可以被视为一个函数$f(x) = Ax$,该函数将$x \in \mathbb R ^k$映射到$y = Ax \in \mathbb R ^n$。
这类函数具有一个特殊性质:它们是线性的。
如果对于所有的$x, y \in \mathbb R ^k$和所有标量$\alpha, \beta$,一个函数$f \colon \mathbb R ^k \to \mathbb R ^n$满足:
则称该函数为线性函数。
你可以验证,当$b$为零向量时,函数$f(x) = A x + b$满足这个性质,而当$b$非零时则不满足。
事实上,已知
再次回顾方程组{eq}la_se。
如果我们比较{eq}la_se和{eq}la_atx,我们可以看到{eq}la_se现在可以更方便地写成
:label: la_se2
y = Ax
我们面临的问题是确定一个向量$x \in \mathbb R ^k$来解{eq}la_se2,其中$y$和$A$是已知的。
这是一个更一般问题的特例:找到一个$x$使得$y = f(x)$。
给定任意函数$f$和一个$y$,是否总是存在一个$x$使得$y = f(x)$?
如果存在,它是否总是唯一的?
这两个问题的答案都是否定的,如下图所示
---
tags: [output_scroll]
---
def f(x):
return 0.6 * np.cos(4 * x) + 1.4
xmin, xmax = -1, 1
x = np.linspace(xmin, xmax, 160)
y = f(x)
ya, yb = np.min(y), np.max(y)
fig, axes = plt.subplots(2, 1, figsize=(10, 10))
for ax in axes:
# Set the axes through the origin
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
ax.set(ylim=(-0.6, 3.2), xlim=(xmin, xmax),
yticks=(), xticks=())
ax.plot(x, y, 'k-', lw=2, label='$f$')
ax.fill_between(x, ya, yb, facecolor='blue', alpha=0.05)
ax.vlines([0], ya, yb, lw=3, color='blue', label='$f$的值域')
ax.text(0.04, -0.3, '$0$', fontsize=16)
ax = axes[0]
ax.legend(loc='upper right', frameon=False)
ybar = 1.5
ax.plot(x, x * 0 + ybar, 'k--', alpha=0.5)
ax.text(0.05, 0.8 * ybar, '$y$', fontsize=16)
for i, z in enumerate((-0.35, 0.35)):
ax.vlines(z, 0, f(z), linestyle='--', alpha=0.5)
ax.text(z, -0.2, f'$x_{i}$', fontsize=16)
ax = axes[1]
ybar = 2.6
ax.plot(x, x * 0 + ybar, 'k--', alpha=0.5)
ax.text(0.04, 0.91 * ybar, '$y$', fontsize=16)
plt.show()
在第一个图中,由于函数不是一一对应的,存在多个解,而 在第二个图中,由于$y$在函数$f$的值域之外,所以没有解。
我们能否对{eq}la_se2中的$A$施加条件来避免这些问题?
在这种情况下,最重要的是要认识到表达式$Ax$对应于$A$的列的线性组合。
具体来说,如果$a_1, \ldots, a_k$是$A$的列,那么
因此,$f(x) = Ax$的值域恰好是$A$的列的张成空间。
我们希望值域足够大,以包含任意的$y$。
你可能记得,我们希望张成空间足够大的条件是{ref}线性独立<la_li>。
一个令人高兴的事实是,$A$的列的线性独立性也能给我们唯一性。
确实,根据我们{ref}之前的讨论 <la_unique_reps>,如果${a_1, \ldots, a_k}$是线性无关的,且$y = Ax = x_1 a_1 + \cdots + x_k a_k$,那么不存在$z \not= x$满足$y = Az$。
让我们讨论更多细节,首先从$A$是$n \times n$矩阵的情况开始。
这是我们熟悉的方程个数等于未知数个数的情况。
对于任意的$y \in \mathbb R ^n$,我们希望找到唯一的$x \in \mathbb R ^n$使得$y = Ax$。
根据上面的观察,如果$A$的列向量是线性无关的,那么它们的张成空间,也就是函数$f(x) = Ax$的值域,就是整个$\mathbb R ^n$。
因此总是存在$x$使得$y = Ax$。
而且,这个解是唯一的。
特别地,以下陈述是等价的:
-
$A$ 的列向量是线性无关的。 - 对于任意$y \in \mathbb R ^n$,方程$y = Ax$有唯一解。
列向量线性独立的性质有时被表述为具有满秩列。
我们能给出某种形式的解吗?
如果$y$和$A$是标量且$A \not= 0$,那么解为$x = A^{-1} y$。
在矩阵的情况下也有类似的表达式。
特别地,如果方阵$A$具有满秩列,那么它就有一个乘法逆矩阵
因此,如果我们用$A^{-1}$左乘等式$y = Ax$的两边,就得到$x = A^{-1} y$。
这就是我们要找的解。
每个方阵都有一个与之唯一对应的数值,这个数值被称为矩阵的行列式。行列式的具体计算方法可以在这里找到。
如果矩阵$A$的行列式不为零,我们就说$A$是非奇异的。
关于行列式最重要的事实可能是:$A$是非奇异的,当且仅当$A$具有满列秩。
这为我们提供了一个有用的单一数值,用来概括一个方阵是否可逆。
这是$n \times k$的情况,其中$n > k$。
这种情况在许多场合都非常重要,尤其是在线性回归中(其中$n$是观测数量,$k$是解释变量的数量)。
对于任意的$y \in \mathbb R ^n$,我们寻找一个$x \in \mathbb R ^k$使得$y = Ax$。
在这种情况下,解的存在性是极不可能的。
在不失一般性的情况下,让我们着重讨论$A$的列向量线性独立的情况。
由此可知,$A$的列向量张成的空间是$\mathbb R ^n$中的一个$k$维子空间。
这个张成空间不太"可能"包含任意的
要理解原因,请回顾{ref}上图<la_3dvec>,其中
想象一个任意选择的
从某种意义上说,这个可能性一定很小,因为这个平面的"厚度"为零。
因此,在
然而,我们仍然可以寻求最佳近似解,例如,找到一个使距离
要解决这个问题,可以使用微积分或正交投影理论。
解为
这是
在这种情况下,要么没有解,要么有无穷多个解 --- 换句话说,解永远不会是唯一的。
例如,考虑
因此,$A$ 的列由
这组向量永远不可能线性独立,因为可以找到两个向量就能张成
(例如,使用标准基向量)
由此可知,其中一列是其他两列的线性组合。
例如,假设
那么如果
换句话说,解不具有唯一性。
下面展示如何使用SciPy的linalg子模块求解线性方程。
所有这些程序都是经过时间检验和高度优化的FORTRAN代码的Python前端接口
A = ((1, 2), (3, 4))
A = np.array(A)
y = np.ones((2, 1)) # 列向量
det(A) # 检查A是非奇异的,因此是可逆的
A_inv = inv(A) # 计算逆矩阵
A_inv
x = A_inv @ y # 解
A @ x # 应该等于y
solve(A, y) # 产生相同的解
请注意我们可以通过inv(A) @ y或使用solve(A, y)来求解$x = A^{-1} y$。
后一种方法使用了不同的算法(LU分解),该算法在数值上更加稳定,因此几乎总是应优先选用。
如果需要求解最小二乘解$\hat x = (A'A)^{-1}A'y$,可以使用scipy.linalg.lstsq(A, y)。
(la_eigen)=
设$A$是一个$n \times n$的方阵。
如果$\lambda$是一个标量,且$v$是$\mathbb R ^n$中的非零向量,满足
则我们称$\lambda$是$A$的特征值,而$v$是特征向量。
换句话说,特征向量是那些经过线性变换$f(x) = Ax$后,方向保持不变、只是长度发生变化的向量。
下图中,蓝色箭头表示两个特征向量,红色箭头表示它们经过变换后的结果。
可以看到,每个特征向量$v$经过变换后得到的$Av$都与原向量指向相同的方向,只是长度改变了
---
tags: [output_scroll]
---
A = ((1, 2),
(2, 1))
A = np.array(A)
evals, evecs = eig(A)
evecs = evecs[:, 0], evecs[:, 1]
fig, ax = plt.subplots(figsize=(10, 8))
# 设置通过原点的坐标轴
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
ax.grid(alpha=0.4)
xmin, xmax = -3, 3
ymin, ymax = -3, 3
ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax))
# 绘制每个特征向量
for v in evecs:
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(facecolor='blue',
shrink=0,
alpha=0.6,
width=0.5))
# 绘制每个特征向量的像
for v in evecs:
v = A @ v
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(facecolor='red',
shrink=0,
alpha=0.6,
width=0.5))
# 绘制它们所在的直线
x = np.linspace(xmin, xmax, 3)
for v in evecs:
a = v[1] / v[0]
ax.plot(x, a * x, 'b-', lw=0.4)
plt.show()
特征值方程等价于
这反过来又等价于行列式为零。
因此,要找到所有特征值,我们可以寻找使
这个问题可以表示为求解一个关于
这进而意味着在复平面上存在
关于方阵
-
$A$ 的行列式等于所有特征值的乘积。 -
$A$ 的迹(主对角线上元素的和)等于所有特征值的和。 - 如果
$A$ 是对称矩阵,那么它的所有特征值都是实数。 - 如果$A$是可逆的,且$\lambda_1, \ldots, \lambda_n$是它的特征值,那么$A^{-1}$的特征值是$1/\lambda_1, \ldots, 1/\lambda_n$。
第一个陈述的一个推论是:矩阵可逆当且仅当它的所有特征值都不为零。
使用SciPy,我们可以按如下方式求解矩阵的特征值和特征向量
A = ((1, 2),
(2, 1))
A = np.array(A)
evals, evecs = eig(A)
evals
evecs
注意 evecs 的列是特征向量。
由于特征向量的任意标量倍数都是具有相同特征值的特征向量(可以验证),eig 程序会将每个特征向量的长度归一化为1。
有时考虑广义特征值问题很有用,即对给定的矩阵
这可以通过 SciPy 中的 scipy.linalg.eig(A, B) 求解。
当然,如果
我们通过简要介绍几个其他重要概念来结束讨论。
回顾几何级数的常用求和公式,即
如果
这个思想在矩阵环境中有一个推广形式。
(la_mn)=
设
右边的范数是普通的向量范数,而左边的范数是矩阵范数——在这种情况下,称为谱范数。
例如,对于方阵
(la_neumann)=
设
换句话说,$A^k$ 是
诺伊曼定理指出:如果对某个
:label: la_neumann
(I - A)^{-1} = \sum_{k=0}^{\infty} A^k
(la_neumann_remarks)=
Gelfand公式告诉我们,对于任意方阵
这里
作为Gelfand公式的结果,如果所有特征值的模都严格小于1,
则存在一个
在这种情况下,{eq}la_neumann 是有效的。
设
我们称
- 如果对于每个
$x \in \mathbb R ^n \setminus {0}$ ,都有$x' A x > 0$ ,则称矩阵为正定的 - 如果对于每个
$x \in \mathbb R ^n$ ,都有$x' A x \geq 0$ ,则称矩阵为半正定的或非负定的
负定和半负定矩阵有类似的定义。
值得注意的是,如果
(la_mcalc)=
以下公式在许多经济学场景中都很有用。设
-
$z, x$ 和$a$ 都是$n \times 1$ 向量 -
$A$ 是$n \times n$ 矩阵 -
$B$ 是$m \times n$ 矩阵,$y$ 是$m \times 1$ 向量
则
$\frac{\partial a' x}{\partial x} = a$ $\frac{\partial A x}{\partial x} = A'$ $\frac{\partial x'A x}{\partial x} = (A + A') x$ $\frac{\partial y'B z}{\partial y} = B z$ $\frac{\partial y'B z}{\partial B} = y z'$
下面的{ref}la_ex1要求你应用这些公式。
scipy.linalg子模块的文档可以在这里找到。
计量经济学理论的第2章和第3章包含了与上述内容类似的线性代数讨论,并附有已解答的练习题。
如果你不介意稍微抽象的方法,{cite}Janich1994是一本不错的中级线性代数教材。
:label: la_ex1
设$x$是一个给定的$n \times 1$向量,考虑以下问题
受限于线性约束
其中
-
$P$ 是一个$n \times n$矩阵,$Q$是一个$m \times m$矩阵 -
$A$ 是一个$n \times n$矩阵,$B$是一个$n \times m$矩阵 -
$P$ 和$Q$都是对称且半正定的
(
解决这个问题的一种方法是构建拉格朗日函数
其中
尝试应用上面给出的关于二次型和线性型求导的公式,得到关于
证明这些条件意味着:
-
$\lambda = - 2 P y$ 。 -
$u$ 的最优选择满足$u = - (Q + B' P B)^{-1} B' P A x$ 。 - 函数
$v$ 满足$v(x) = - x' \tilde P x$ ,其中$\tilde P = A' P A - A'P B (Q + B'P B)^{-1} B' P A$ 。
正如我们将看到的,在经济学背景下,拉格朗日乘子通常是影子价格。
如果我们不关心拉格朗日乘数,我们可以将约束条件代入目标函数,然后仅对$u$最大化$-(Ax + Bu)'P (Ax + Bu) - u' Q u$。你可以验证这会得到相同的最大值。
:class: dropdown
我们有一个优化问题:
满足约束条件:
其中基本条件为:
-
$P$ 是一个对称且半正定的$n \times n$矩阵 -
$Q$ 是一个对称且半正定的$m \times m$矩阵 -
$A$ 是一个$n \times n$矩阵 -
$B$ 是一个$n \times m$矩阵
相关的拉格朗日函数是:
第1步:
对拉格朗日方程关于y求导并令其导数等于零得到:
因为P是对称的。
因此,关于y的拉格朗日方程最大化的一阶条件意味着
第2步.
对拉格朗日方程关于u求导并令其导数等于零得到
代入$\lambda = -2 P y$得到
将线性约束$y = Ax + Bu$代入上式得到
这是关于u的拉格朗日方程最大化的一阶条件。
因此,u的最优选择必须满足
这是由拉格朗日方程的一阶条件定义得出的。
第3步.
将约束代入目标函数,重写我们的问题,得到
由于我们知道u的最优选择满足$u = -(Q + B'PB)^{-1}B'PAx$,那么
$$
v(x) = -(Ax+ B u)'P(Ax+B u) - u'Q u ,,,, 其中 ,,,, u = -(Q + B'PB)^{-1}B'PAx $$
计算函数
为简化起见,令
对于第二项
注意到项
对于第三项
因此,第二项和第三项的和为
这意味着
因此,优化问题的解