PYTHON替代MATLAB在线性代数学习中的应用(使用Python辅助MIT 18.06 Linear Algebra学习)

前言

MATLAB一向是理工科学生的必备神器，但随着中美贸易冲突的一再升级，禁售与禁用的阴云也持续笼罩在高等学院的头顶。也许我们都应当考虑更多的途径，来辅助我们的学习和研究工作。 虽然PYTHON和众多模块也属于美国技术的范围，但开源软件的自由度毕竟不是商业软件可比拟的。

本文是一篇入门性文章，以麻省理工学院（MIT) 18.06版本线性代数课程为例，按照学习顺序介绍PYTHON在代数运算中的基本应用。 介绍PYTHON代数计算的文章非常多，但通常都是按照模块作为划分顺序，在实际应用中仍然有较多困扰。 而按照代数课程学习的顺序，循序渐进，集注在最常用和最实用的功能上，比较符合典型的应用逻辑。可以用较低的门槛逐步完成PYTHON在线性代数中各项功能的学习和应用。 MIT 2020版本的线性代数课程也已发布，但基本是在18.06版本上的修正。Gilbert教授的年龄已经很大，只录制了一个5节课的串讲。所以系统性还是18.06版本更为完整。 很讽刺是吧，课程本身也是美国的-_-#。阿Q一下吧，就当是“师夷长技以制夷”。

首先给出几个相关链接： MIT 18.06 Linear Algebra课程主页 B站完整版34讲Gilbert教授课程视频 配套第三版线性代数教材(百度网盘) 提取码：uhvc 最新发行的教材是第5版，建议听课时使用配套的第3版教材。课程完成后，把第5版教材作为辅助读物。不然在章节、内容方面会碰到很多困惑。

版本选择

PYTHON版本的选择现在已经没有什么困惑了，PYTHON2停止了支持，PYTHON3现在是必选项。我是用Mac电脑，通常使用brew安装PYTHON3，每次有新版本的时候执行brew upgrade会自动升级。不使用内置的PYTHON3是为了防止安装很多扩展库的时候会有系统完整性检查导致的不兼容，不过只是跑一跑数学运算的话倒也不用担心这些。 Linux各发行版各发行版是Python的开发环境，所以内置的PYTHON3就很好用。使用apt/yum等包管理工具升级的时候会自动完成版本维护。 PYTHON在Windows/Linux/Mac等各平台上兼容性非常好，特别是在数学计算方面基本不用担心互相之间的通用问题。

计算模块方面则有了很多的选择，常见的有NumPy/SciPy/SymPy。 其中在数值计算方面，NumPy应用非常广泛，特别是TensorFlow/PyTorch等机器学习平台也把NumPy当做默认支持之后。所以本文会把NumPy当做一个选择。 在课程学习和理论研究方面，符号计算更为重要。SymPy在这方面有比较多的优势，所以本文会把SymPy当做另外一个选择。 SciPy以及还有一些小众计算模块同样非常优秀，但限于篇幅，本文中只好做一些取舍。

在PYTHON3环境下安装NumPy/SymPy模块的方法很简单：

如果碰到麻烦，一般是因为网络速度造成的。特别是默认的国外软件源。修改软件源到国内的服务器会提高不少下载速度，方法是修改文件~/.pip/pip.conf，默认是没有这个文件的，要自己建立~/.pip目录和新建对应的文本文件，内容为：

这里使用了清华大学的镜像服务器。 以上是在Linux/Mac之上的操作方法。Windows用户，虽然PYTHON3本身没有兼容问题，但还是建议你使用Windows10内置的Linux子系统来学习。因为Mac/Linux下Python的资料更为丰富，能让你节省很多时间。

矩阵的表达

在Pyhton中使用扩展库，首先要做引用，比如引入NumPy库：

意思是引用numpy计算库，并重新命名为np。使用numpy中的方法时，首先要以“np.”开头。 SymPy库的引用，通常会直接从中将所有资源直接引用到当前作用域，像使用原生方法一样使用SymPy中定义的方法，这也是SymPy官方推荐的：

出于个人习惯，我还是更喜欢同使用NumPy一样使用SymPy:

虽然因此所有的SymPy的方法都要冠以“sp.”前缀，但这样不容易造成混淆从而导致错误。

以下内容大致对应课程（MIT 18.06 Linear Algebra课程，以下简称课程）第一、二讲的内容。 在线性代数中，主要涉及3种数据类型，常量、标量（Scalar)、向量(Vector)、矩阵（Matrix)。 无论NumPy还是SymPy，都直接使用了基本Python类型作为标量，比如：

而对于向量和矩阵，处理方法则有很大区别，下面先讲述NumPy中的方法。

假设我们有v1、v2两个向量，及A、B两个矩阵：

\[v1 = \left[\begin{matrix}1\\2\\\end{matrix}\right] \]

\[v2 = \left[\begin{matrix}3\\4\\\end{matrix}\right] \]

\[A = \left[\begin{matrix}1&2\\3&4\\\end{matrix}\right] \]

\[B = \left[\begin{matrix}5&6\\7&8\\\end{matrix}\right] \]

1. 首先，NumPy接受Python原生的数组当做向量和矩阵

1. 其次，NumPy内置的数组类型（Array)也可以用来表示向量和矩阵

1. 正式的矩阵表示方法，使用NumPy内置的Matrix类型

第1、2种方法，虽然在很多情况下都能正常的使用，但都算不上规范化的矩阵表示方法。特别是对于向量的表示，向量本来是纵向的，代表矩阵中的一列。但在方法1、2中，都成为了横向的。这很容易造成概念的混淆，和计算中的错误。 当然Python内置的列表类型以及NumPy内置的列表类型并非不能使用，实际上它们在计算速度上有非常明显的优势。简单说就是功能少的类型，往往有高的速度。所以渐渐的我们能看到很多需要速度的运算，还是会使用np.array类型操作。实际上对各种类型熟悉了之后，有了自己的习惯和原则，什么时候用什么类型，你自然会有自己的标准。 为了让大家对这种差异有更清晰的认识，这里举几个例子，也顺便看一看最基本的矩阵计算：

以上的示例可以明显看出，对于Python内置数组类型，并没有定义对应的矩阵操作，所以不能直接用于线性代数的计算。NumPy的很多方法都接受使用Python内部数组作为参数来表达向量和矩阵，所以给人的印象，这些类型之间没有什么区别。 NumPy内置的数组类型和矩阵类型，在简单运算中都能得到正确的结果，可以用于常用的计算。但实际上很多高级函数及算法，对两种类型的处理仍然存在很大区别，就类似示例中出现的矩阵乘法。所以在彻底了解之前，不建议使用np.array类型当做矩阵类型来使用。否则在复杂的项目中，很多莫名其妙的计算错误会让你排错工作异常复杂。

NumPy还提供了一种更方便的方法来定义向量和矩阵，是我当前最常用的方法：

熟悉MATLAB的同学应当开心了，这跟MATLAB中定义矩阵的方法完全一样，算是Python环境中最方便的方法。

在线性代数课程中，经常会需要选取一个典型矩阵，做一些计算的练习。课堂上Gilbert教授经常随手就可以举出一个矩阵的例子，并且各行列线性无关，而我们往往很难做到这一点。这时候可以使用随机数初始化矩阵的方法：

与此类似的，还有初始化一个全0矩阵、全1矩阵、单位方阵I的方法，如果你打算用程序逻辑建立一个矩阵，那这些往往是开始的第一步：

下面看看SymPy定义向量、矩阵的方法。

作为符号计算的代表，SymPy的计算结果通常都是公式形式，所以SymPy专门提供了LaTeX的输出方式：

这种输出格式对通常的程序没有什么意义。但如果是用于论文写作的话，可以直接拷贝到LaTex编辑器，成为一个精致的公式输出。就类似本文中的公式，通常也是采用LaTeX格式输入。

求解线性方程

这也是课程第一、二讲中的内容。方程组是矩阵的起源，也是矩阵最初的目的。以课程中的方程组为例：

\[\left\{ \begin{array}{l} 2x-y=0\\ -x+2y=3 \end{array} \right. \]

可以得到矩阵A及向量b:

\[A = \left[\begin{matrix}2&-1\\-1&2\\\end{matrix}\right] \]

\[b = \left[\begin{matrix}0\\3\\\end{matrix}\right] \]

\[Ax=b \]

这里的x实际代表两个未知数组成的向量：

\[x = \left[\begin{matrix}x\\y\\\end{matrix}\right] \]

使用NumPy解方程组示例：

使用SymPy解方程组示例：

作为符号计算的优势，SymPy中可以定义未知数符号之后，再使用跟NumPy中同名的方法solve()来直接对一个方程组求解，但那个不属于本文的主题范畴，所以不做介绍。有兴趣的话也可以参考这篇老博文《从零开始学习PYTHON3讲义（十一）》。 SymPy跟NumPy语法差异还是比较大的，使用中需要特别注意。两个软件包，虽然都是Python中的实现，但并不是由同一支开发团队完成的。所以这种差异感始终是存在的。

矩阵乘法和逆矩阵

这是课程第三讲的内容，其中矩阵同矩阵的乘法运算在前面一开始就演示过了，对于手工计算来讲，这是最繁琐的部分。而对于Python,这是最简单的部分。 矩阵的逆在线性代数中会频繁出现，经常用到，两个软件包中都已经有了内置的方法。

下面在同一个代码块中分别演示NumPy和SymPy的方法：

上面代码非常明显的体现出了NumPy数值计算和SymPy符号计算的不同。前者会因精度所限有极小的误差，而后者通常能保持美观的整数数字。但前者的数字可以直接应用到机器学习等业务系统。而后者是对人的理解更有益，归根结底还是符号，不能当做数值使用。 好在Python之中，如果不考虑转换速度，不同模块之间共享数据非常容易。前面的演示中已经有了将NumPy矩阵转换为SymPy矩阵，以及将SymPy的计算结果转换到NumPy的实例。这对用户来说，是非常方便的。

矩阵的LU分解

课程第四讲重点讲解了矩阵的LU分解。对于一个给定矩阵A,可以表现为一个下三角矩阵和一个上三角矩阵乘积的形式：

\[A=LU \]

其中上三角矩阵U是求解方程组的初步中间产物。由这一步开始，逐步求解靠后的主元，再回代至方程，以求解更多的未知数主元。重复这个步骤，直到完成所有未知数的求解。 NumPy中，并没有提供矩阵的LU分解功能。可能是因为觉得L、U矩阵用途并不是那么广泛，并且可以直接用方程求解来替代。 如果需要用到的话，通常方式是使用其它软件包替代，比如SciPy。 这里也提供一个架构于NumPy之上的子程序，来完成LU分解的功能。子程序内部是将矩阵类型转换为数组类型，从而方便遍历。接着是使用手工消元相同的方式循环完成LU分解。 需要说明的是，这类附带了子程序的Python片段，建议还是保存到一个文本文件中，以脚本方式执行。在交互式方式下很容易出现各种错误。

程序执行可以获得类似这样的输出：

偏重于计算分析的SymPy则直接内置了LU分解功能，对速度不敏感的场合，使用SymPy做LU分解，再转换到NumPy矩阵也是一个好办法：

LU分解那一行的代码，使用下划线忽略的部分是函数返回的行交换矩阵。 在消元过程中，对应主元位置如果为0的话会导致消元失败，此时会产生行交换。这种情况下，会由单位矩阵I变换而来的行交换矩阵先同矩阵A相乘，从而将主元为0的行交换到其它行，保证消元的顺利进行。 使用Python辅助解方程，这些步骤都是很少需要手工操作了，如果有必要，就自行赋值给矩阵变量保留吧。 顺便提一句，讲到置换矩阵的时候，教授还提到了对于一个n*n的方阵，置换矩阵可能有多少种呢？答案是n!，也就是n的阶乘。 在Python内置的数学库、NumPy、SymPy中，都有求阶乘的函数：

第四讲还介绍了矩阵的转置，这是线性代数中使用极为高频的功能：

简化行阶梯矩阵、0空间基、特解、通解

课程第五至第十讲围绕着矩阵的四个基本空间展开。推导和计算很多，但都是基础线性组合，用Python当成计算器就够用了。 在空间维度判断方面，我们倒是能帮上一些小忙，就是计算矩阵的轶。 矩阵的行空间、列空间轶都是相同的。0空间维度是n-r,左0空间维度是m-r。

如果方程组满轶，也就是方程组有解的情况下，开始一节介绍的解线性方程组很不错。 非满轶的情况，求方程组的特解和通解。将矩阵化简为“简化行阶梯矩阵（Reduced Row Echelon Form）”会非常有用。可惜的是，NumPy依然没有提供内置的支持。自己实现的话，代码还挺长的，远不如使用现有的软件包更方便。所以在这里我们主要看SymPy的实现：

函数返回一个RREF矩阵还有一个元组，元组指示RREF矩阵中主元所在的列。这个元组是非常必要的，在第二个例子中就能明显看出，主列并不一定是从左到右相邻排列的。 此时，可以通过RREF最下面的全0行跟方程组b向量的情况判断函数可解性。以及根据自由变量F子矩阵的情况获得方程的0空间解。 当然，如同前面的解方程一样，SymPy中直接提供了函数获取0空间解。RREF这些中间过程主要用于分析学习，对于使用工具自动解方程意义并不大：

方程组的通解包括特解和0空间基两部分。前面获得的是0空间的基。特解则需要方程式右侧向量的配合：

Bs矩阵同b向量的组合获得一个有限集的解，那么这个解中的tau0/tau1是什么意思呢？ 参考前面的rank计算或者rref矩阵，我们知道Bs矩阵有两个自由变量（由n-r得来），tau0/tau1就是这两个自由变量。这也是因为我们没有定义未知数符号所导致的自动命名。如果需要，我们可以定义x1/x2...这样的未知数。不过这不是我们的重点，请忽略这个命名。 方程的特解是当自由变量为0的时候，方程的解。因此将tau0/tau1都设为0，化简一下，很容易得到方程的特解为：      (-2,0,3/2,0)。 再结合上面计算的Bs矩阵在0空间的2个基，就是方程组的通解:

\[X_{Complete}= \left[\begin{matrix}-2\\0\\\frac32\\0\\\end{matrix}\right]+ k1\left[\begin{matrix}-2\\1\\0\\0\\\end{matrix}\right]+ k2\left[\begin{matrix}2\\0\\-2\\1\\\end{matrix}\right] \]

点积、获取指定行向量和列向量、正交判定

点积也称作点乘、内积，是向量、矩阵中最常用的一种运算。NumPy和SymPy中都有支持。

矩阵运算中，使用一个向量的转至乘另外一个向量，或者乘自己用于求平方，都是非常常用的使用方法。在使用NumPy做运算的时候要特别注意一点，这样点积的结果仍然是一个矩阵，只是1维*1维。 在线性代数课程上，都会直接把这个点积结果继续用于计算，但在使用NumPy的时候，要特别注意应当将其转换为浮点数，然后再用于计算。不然会出现矩阵维度不符的错误。示例如下：

SymPy作为主要用于实验分析的符号计算工具，点积运算的结果直接就是可继续用于计算的数字，不需要另行转换。

获取矩阵的特定行向量和列向量，在NumPy/SymPy中都是重载了Python语言的列表（数组）操作符，所以方法都是相同的。 需要注意的是在数学中，矩阵行列的计数通常从1开始，第1行、第2行...第1列、第2列。而在Python中，遵循了计算机语言一贯的习俗，是从0开始计数的。Python中矩阵的第0行，就相当于通常数学课程上所说的第1行。 先来看获取矩阵中特定元素的方法：

获取行向量、列向量，相当于获取矩阵某一行或者某一列所有的数据。在Python中，使用':'字符放置在行、列参数的位置，就代表获取完整行或者列的数据：

点积和正交判断是在课程第十四讲中引入的。 判断两个向量是否正交，就是用一个向量的转置，点积另外一个向量。相互正交的向量，点积结果为0。上面的例子说明，我们随意定义的矩阵，前两列并不正交。 单位矩阵I的每一行、每一列都是正交的，我们测试一下：

此外在NumPy和SymPy的运算符重载中，乘法运算符'*'直接就定义为了点积运算，是可以直接使用的：

方程组的最优解

内容同样来自课程第十四讲。 在实际的应用中，方程组的数据来源经常是测量的结果。在一组实验中，测到了多组结果，这代表方程有多行。但因为测量误差以及干扰因素，这些不准确的测量值所形成的方程组，往往因为误差导致的矛盾因素是无解的。 这时候，通过计算测量数据到方程组矩阵列空间的投影信息，形成新的方程组，可以得到最接近真实结果的答案，这就是最优解。 对于一个原始方程：

\[Ax=b \]

其最优解方程为：

\[A^TA\hat{x}=A^Tb \]

求得的\(\hat{x}\)就是方程的最优解。它并不是原来的x,而是最接近合理值的近似解，所以称为最优解。 下面使用SymPy为例演示方程组求解最优解,NumPy可以使用同样的方法：

投影矩阵

投影矩阵的概念来自课程第十五讲。 使用投影矩阵公式可以求得矩阵A的投影矩阵：

\[P=A(A^TA)^{-1}A^T \]

下面以NumPy为例，演示计算投影矩阵：

正交矩阵和正交化法

这部分内容来自课程第十七讲。 按照教授的说法，标准正交矩阵是能得到的最好的矩阵，有很多优良的性质，便于计算和分析。 标准正交矩阵每一列都互相垂直，模长为1。通常把标准正交矩阵记为Q。 但很可惜，通常的矩阵都不是标准正交矩阵。课程中介绍了格拉姆-施密特(Graham-Schmidt)正交化法，将一个列满轶的矩阵A，转换为一个由标准正交向量组构成的矩阵Q。 SymPy内置了这个算法，用于将一组线性无关的向量正交化，来看看示例：

这个小程序段需要单独保存为一个脚本来执行，输出因为SymPy符号计算的特点，会变得极为复杂。这种复杂主要来自于标准化除以模长所导致的分数化。

毕竟不是编程的课程，所以虽然是很短一个小程序，非IT专业的同学看起来可能也会觉得晕。这是由于SymPy中内置的格拉姆-施密特算法主要用于处理向量所导致的。我们不得不把矩阵变为向量，完成正交化后，再转换回矩阵。 实际上有更好的办法，就是使用QR分解。QR分解计算起来更麻烦，在课程中并没有介绍，不过还是老话，计算机最不怕的就是清晰的计算。 QR分解的大意是，任何一个列满轶的矩阵A,都可以分解为一个标准正交向量Q和一个上三角矩阵R的乘积形式。上三角矩阵前面见过，就是我们使用高斯消元的中间步骤产物U。 SymPy和NumPy中都内置了QR分解算法，请看示例：

行列式、伴随矩阵、特征值、特征向量

这几个概念可以说是线性代数的核心，因为计算太复杂，贯穿了多讲内容，从第十八讲一直延续到了第二十一讲。 其中为了降低行列式的计算量，还穿插了代数余子式。但计算机的发展让这些复杂计算都成为了一行函数的事情，所以很多基本的加法、乘法的运算，我们就忽略掉了。 这部分没有太多可说的，直接用示例来说明吧：

对角矩阵和对角化

这部分内容来自课程第二十二讲。 矩阵A对角化为Λ的公式为：

\[Λ = S^{-1}AS \]

S是矩阵A特征向量所组成的矩阵。 下面用numpy示例一下：

嗯，为了验证课程中的公式，故意搞复杂了点。这样的计算其实完全没有必要，对角化矩阵实际就是矩阵特征值排列在对角线所组成的矩阵。所以实际上把特征值乘单位矩阵I，转化到对角线就好了：

SymPy也可以使用对角化公式计算，但SymPy计算的特征向量需要自己解析、组合成矩阵S，有点麻烦。 幸运的是，SymPy直接提供了矩阵对角化的函数：

SymPy的计算结果看上去总是那么工整。既然有了S和Λ，是不是想用SymPy算回到矩阵A验证一下？不不不，你不会想那么做的，不信我给你练练：

看起来很晕是吧？ 符号计算有符号计算的优点，但缺点也那么显而易见。速度慢就不说了，复杂计算的时候，表达式化简能力，特别是灵活程度毕竟不能同人相比，这就是一个很典型的例子。这样的结果，肯定是还不如用NumPy计算的近似值。 怀疑计算出了错？那倒肯定没有，我们把上面矩阵用NumPy计算出最终数值看一下：

跟最初的矩阵A是吻合的，毋庸置疑。

矩阵乘幂、矩阵法求斐波那契数列、绘图

同样来自课程第二十二讲，对角化矩阵的一种典型应用就是简化矩阵的幂运算。 对于一个高维矩阵的高次幂来讲，如果手工计算，其复杂程度是可想而知的。而通过对角化后的矩阵，矩阵幂的运算可以简化很多：

\[A^k = SΛ^kS^{-1} \]

使用计算机之后，这种简化手段意义显得不再那么显著。但这种思想还是非常有帮助。 比如对于计算序列值的时候，这成为了一种全新的思路，类似符合\(u_{k+1} = Au_k\)这样公式的数字序列，都可以用这个思路来计算。在\(u_0\)已知的情况下，公式可以变形为\(u_k=A^ku_0\)。 以电脑编程入门必学的斐波那契数列生成为例，通常是使用循环顺序生成（此处略去程序）。 使用矩阵运算的思想，联立方程（方程请参考原课程）可以得到矩阵：

\[A = \left[\begin{matrix}1&1\\1&0\\\end{matrix}\right] \]

以及初始向量为：

\[u_0 = \left[\begin{matrix}1\\1\\\end{matrix}\right] \]

想得到序列中任意一个数值，无需循环，一次计算就可以直接得到。来看一下获取斐波那契数列的代码片段：

把上面代码保存为脚本文件，执行后获得输出为：

线性代数是研究向量和空间的学科，绘图能够在很大程度上帮助问题的思考。前面一直没有合适的例子，在这里牵强的引入一下，就是将计算结果用绘图的方式展示出来。 接着上面程序代码继续在后面加入：

程序执行前首先要安装绘图软件包pip3 install matplotlib。安装完成后再次执行程序，除了得到上面相同的4行输出外，还会绘出如下的曲线，表现了斐波那契数列序号同数值之间的关系：

在Win10的Linux子系统使用绘图功能的话，还需要配置X11协议的远程显示，请参考这篇文章：win10配置linux子系统使用python绘图并显示

马尔科夫矩阵

内容来自课程第二十四讲。 马尔科夫矩阵也叫状态转移矩阵，用于表现在总数不变的情况下，随着时间的延续，不同节点之间某项指标的动态情况。 课程中的例子是讲伴随着硅谷的崛起，加州人口逐渐增加。而相邻麻省人口受加州经济吸引，移居至加州，从而导致本省人口减少。这种情况持续一段时间之后，最终形成一个稳态。例子中的数字显然就是示意性。 马尔科夫矩阵代表了使用矩阵思想，在具体事例中的应用。 下面代码使用课程中的数据做了一个计算推演，并且绘制了曲线图供参考。代码需要保存为脚本文件执行：

程序执行绘制的曲线如下：

程序在绘图中，使用了中文的图例注释。如果不做任何处理，出现的会是乱码。plt.rcParams['font.sans-serif']=['Songti SC']这一行代码，就是将默认绘图字体名sans-serif指定到使用系统宋体，从而正常显示中文。在不同的电脑上，要根据自己的电脑字体名称设置，选择一个替换。

对称矩阵、复矩阵

这部分内容来自课程第二十五、二十六讲。 对于实数矩阵来说，对称矩阵就是转置与自身相同的矩阵，判定起来很容易。 实对称矩阵还有特征值全部为实数、特征向量相互垂直两个重要特性。

SymPy的相关操作前面都已经出现过，此处不再重复。

复矩阵就是元素中存在复数的矩阵。关键是复数如何表达，NumPy中延续了Python中对复数的定义方式；SymPy中定义了自己的虚数符号类。两种方式都离我们日常数学中的习惯区别很大。

对称复矩阵(埃尔米特矩阵)的定义跟实数矩阵有所区别，在复矩阵中，对称是指矩阵做完共轭、转置操作后，同本身相等。这也意味着，在对称复矩阵对角线上的元素必须都是实数。否则不可能做到共轭后与自身相同。 复矩阵组成的正交矩阵称为酉矩阵。

正定矩阵

正定矩阵是课程第二十七讲的内容。 首先经常用到的是正定矩阵的判定。课程中教授提了一个问题：

\[A = \left[\begin{matrix}2&6\\6&c\\\end{matrix}\right] \]

A矩阵定义如上，右下角c取值是多少，使得A矩阵成为正定矩阵。 老师给了几个人工判定的标准：

1. 矩阵为对称方阵。
2. 所有特征值为正。
3. 所有主元为正。
4. 从左上角开始的子对称矩阵行列式为正。
5. 对于任意非零向量x,xᵀAx的结果为正。

对于SymPy来讲比较容易，内置提供了正定矩阵判定的方法。NumPy没有内置此种功能，但可以根据上面的标准，用一小段程序来判断，难度也不大。不过NumPy还有一个取巧的办法，NumPy中有矩阵的霍尔斯基分解函数，霍尔斯基分解是要求矩阵为正定矩阵的。如果提供的矩阵参数不是正定矩阵，函数会报错。此时截获这个错误，也可以准确、方便的判定矩阵的正定性。 （霍尔斯基分解不属于本课程内容，这里只是利用它来判断矩阵的正定性，所以不用管其具体功能。） 下面用代码来看看具体方法：

接着以上面的矩阵为例，来实际判定测试一下：

回到老师出的问题，课程中讲过，通过xᵀAx>0来判断矩阵的正定性是最全面、可靠的。 以题中的两维矩阵为例，向量也就是两维，假设为：[x1;x2]，把矩阵A代入公式、展开：

\[x^TAx = 2x_1^2+12x_1x_2+cx_2^2 \]

可以看出，第一个系数2，就是矩阵A左上角的值；第二个系数12是A第1行第2列及第2行第1列的和；第三个系数就是c了。实际上这来自于行列式的计算。由此可判断c>18可以保证矩阵A是正定矩阵。 对于课程的内容我没有什么要补充的。但是在Python的帮助下，如果将上面公式图示出来，肯定可以帮助我们更深入的理解矩阵A中c取值对于矩阵正定性的影响。 因为上面公式有x1/x2两个变量，加上最终整体公式的取值算作一个维度，我们需要绘制的是三维图。

下面程序中，我们分别使用c=7以及c=20，绘制两幅三维图片。程序使用了matplotlib绘图软件包的mpl_toolkits三维绘图模块。这个模块是matplotlib新版本才引入的，所以如果找不到这个模块的话，请升级matplotlib模块。

绘制结果请看图：

在第一张图片中，可以看到当c取值7时，有相当一部分图都位于0（Z轴)之下。 而第二张图片中，c取值20，所有曲线都会在0之上了，代表xᵀAx>0，矩阵是正定矩阵。 绘制的三维图片，可以使用鼠标拖动，从各个角度观察。还可以旋转、缩放、保存为图片文件。Python实在是数学学习和研究的利器。

奇异值分解

这是课程第二十九讲的内容。奇异值分解的公式如下：

\[A = U∑V^T \]

其中，U是AAᵀ矩阵的特征向量形成的标准正交矩阵；V是AᵀA矩阵的特征向量形成的标准正交矩阵；∑则是两个矩阵特征值开根号后形成的对角矩阵。 SVD分解几乎串起了整个线性代数课程的知识点，手工计算的话，还是相当的麻烦。 NumPy中已经内置了奇异值分解的函数：

这次终于碰到了SymPy没有内置对应功能的情况。实际你也看出来了，SVD是典型的数值计算应用，对于符号运算分析的话作用并不大。而且因为运算步骤多，符号计算保留过多的符号操作很容易造成计算溢出，可读性更是没有了保障。 所以在SymPy的官方推荐中，也是使用mpmath运算包完成SVD分解。在新版本的SymPy中，这个包已经分离并且需要单独安装，所以你还不如直接使用NumPy计算了。 上面的计算中，变量s代表了SVD分解之后的∑对角矩阵，实际是AAᵀ矩阵或者AᵀA矩阵特征值再开方的值。使用NumPy做完SVD分解后，直接保存为列表类型。对角阵除了对角线的数据，其它元素都是0，这样保存非常合理。 把列表还原回对角阵，前面介绍了乘单位矩阵I的方法，这里再介绍一个NumPy内置的diag函数，用起来更方便：

最后

本文基本涵盖了MIT18.06版线性代数课程中所需要用到的计算方法。很多章节中，计算量虽然不小，但已经在前面讲过的，就被省略掉了。 希望能为学习线性代数的同学，或者希望从MATLAB迁移至Python的同学带来帮助。 错误疏漏因水平所限难以避免，欢迎指正。