多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

MATLAB线性代数实战:从矩阵基础到方程组求解与矩阵分解

MATLAB线性代数实战:从矩阵基础到方程组求解与矩阵分解 MATLAB 线性代数是 MATLAB 使用频率最高、也最容易被低估的部分。很多人把 MATLAB 当成一个高级计算器看到-*就直接写却忽略了 MATLAB 从诞生之初就以矩阵为最小计算单位。真正进入工程场景后线性方程组的求解方式、矩阵分解的选择、维度不匹配的报错都会成为影响调试效率的关键因素。这份免费 MATLAB 教程的线性代数分篇会从矩阵的基本操作讲起接着覆盖线性方程组求解、矩阵分解、特征值与奇异值再落到最小二乘拟合、图像处理、Simulink 和光学仿真等常见工程场景。学完之后你可以建立一条清晰的 MATLAB 线性代数操作框架知道什么情况用A\b什么情况用eig什么时候该检查条件数什么时候要考虑稀疏矩阵。1. 矩阵在 MATLAB 中不只是数据结构更是计算单位1.1 为什么 MATLAB 以矩阵为最小计算单位MATLAB 的名称来源于 Matrix Laboratory这个背景决定了它的大部分计算逻辑都围绕矩阵展开。即使是一个普通的数值5在 MATLAB 中也可以看作一个1×1的矩阵一个向量是一个1×n或n×1的矩阵一张灰度图像是一个m×n的数值矩阵一组时间序列数据也可以按一列或多列组织成矩阵。这种设计带来一个直接好处运算规则高度统一。线性代数中的矩阵乘法、转置、求逆、分解、特征值计算在 MATLAB 中都有对应的原生函数不需要自己再去写循环。对于工程人员来说理解这一点非常重要因为很多看似复杂的报错最后都能归结为“矩阵维度不匹配”或“把矩阵运算当成了标量运算”。实际项目里最常见的误区是把线性代数当成数学课的独立章节学完就忘。但在 MATLAB 中线性代数不是孤立知识它服务于控制系统、信号处理、图像处理、机器学习、光学仿真、Simulink 建模等场景。只有把矩阵运算作为底层语言才能在看懂别人代码时更快抓住核心。1.2 学习环境与版本选择能跑起来比追求最新版本更重要学习 MATLAB 线性代数环境准备不需要太复杂。MATLAB 是商业软件正式使用需要授权。常见做法是使用学校或公司的授权版本也可以到官网申请试用或者使用 MATLAB Online 这类在线方式。具体安装和授权流程请以 MathWorks 官方信息为准不要盲目相信网上所谓的“最新版安装包”。版本选择上线性代数的基础语法在近几个版本里变化不大。如果你的版本是 R2021a、R2022b、R2023b或者更新的版本本教程中的代码基本都能运行。真正影响运行的是某个具体函数是否属于已安装的工具箱以及你的 MATLAB 是否配置好了编译器或外部接口。安装完成后建议先在命令行窗口确认环境version verversion显示 MATLAB 版本号ver列出已安装的工具箱列表。如果后面运行eig、svd这些函数时提示找不到函数先回到这一步确认工具箱是否完整。环境准备阶段可以按下面这个清单检查检查项说明是否需要MATLAB 主程序能正常打开命令行窗口必需授权状态能运行license查看授权信息必需工作目录进入一个可写的项目目录建议工具箱列表ver查看是否有需要的工具箱按需示例数据内置示例或自行构造数据建议1.3 从“变量”思维切换到“矩阵”思维很多初学者在 Python 或其他语言里习惯了“变量 for 循环”到了 MATLAB 仍然沿用这个习惯。线性代数场景下更推荐的做法是先设计矩阵结构再调用矩阵函数。理解这一点可以先从检查维度开始。看下面这段代码x 1:5; size(x)x并不是一个“数组大小不确定”的东西size(x)会返回1 5说明它是一个行向量。有人想生成列向量会写成y 1:5; size(y)这里要特别小心。1:5在 MATLAB 中并不是(1:5)因为单引号优先级会先作用于5而5仍然是5所以结果依然是一个1×5行向量。如果想得到列向量应该写成z (1:5); size(z)这就是典型的“矩阵思维”问题。不要凭直觉写要时刻问自己当前变量是行向量还是列向量矩阵乘法要求的内维是否匹配输出结果会是什么形状。可以用whos查看变量的具体信息whos x y z这个命令会显示每个变量的名称、大小、字节数、类型。当代码出现Matrix dimensions must agree报错时第一时间运行size或whos往往比逐行读代码更快。2. 从建矩阵到取数据先掌握矩阵的基本操作2.1 创建矩阵的常见方法创建矩阵最直接的方式是使用方括号。同一行元素用空格或逗号分隔不同行用分号分隔A [1 2 3; 4 5 6; 7 8 9];这里A是一个3×3矩阵。如果想创建更复杂的矩阵可以使用 MATLAB 内置函数Z zeros(3, 4); % 3 行 4 列全零矩阵 O ones(2, 5); % 2 行 5 列全一矩阵 I eye(4); % 4×4 单位矩阵 R rand(3, 3); % 3×3 均匀分布随机矩阵 N randn(3, 3); % 3×3 正态分布随机矩阵 D diag([1 2 3]); % 对角矩阵创建矩阵之后要观察它的形状。不要只靠“我觉得它应该是这个维度”用size来确认size(A) size(R)如果要把多个矩阵拼接起来行方向用逗号列方向用分号B [A, A]; % 横向拼接列数翻倍 C [A; A]; % 纵向拼接行数翻倍横向拼接要求两个矩阵行数相同纵向拼接要求列数相同。如果不一致会直接报错这也是Dimensions of arrays being concatenated are not consistent的常见来源。2.2 索引、切片与变形MATLAB 的矩阵索引使用圆括号第一个下标是行第二个下标是列。例如A [1 2 3; 4 5 6; 7 8 9]; a12 A(1, 2); % 第 1 行第 2 列结果为 2 firstColumn A(:, 1); secondRow A(2, :); sub A(1:2, 2:3); % 取第 1 到第 2 行、第 2 到第 3 列A(:)会把矩阵按列顺序拉成一个列向量。这是线性索引的体现MATLAB 的矩阵在内存中按列优先存储因此A(5)访问的是第 5 个元素不是第 5 行第 5 列。变形操作中reshape最常用V 1:12; M reshape(V, 3, 4);reshape要求变换前后元素总数一致。1:12共 12 个元素可以变成3×4或4×3但不能变成3×5。复制平铺可以用repmatrowVec [1 2 3]; rep repmat(rowVec, 2, 1);结果为两行[1 2 3]。repmat会显著扩大内存处理大矩阵时不要滥用先确认是否真的需要。转置也是常见操作。实数矩阵转置直接用A没问题但如果矩阵是复数A是共轭转置而A.才是普通转置。这个细节在信号处理和复数矩阵运算中经常成为隐性 bug。2.3 常用矩阵操作函数速查表日常 MATLAB 线性代数代码中下面这些函数出现频率最高函数作用常见注意点size(A)返回矩阵维度返回的可能是多个值用[m,n] size(A)接收numel(A)返回元素总数等价于m*nnorm(v)计算向量或矩阵范数默认是 2 范数可指定norm(A,fro)det(A)计算行列式浮点误差下小数值不等于奇异rank(A)计算数值秩与容差有关不要当作绝对精确trace(A)计算迹方阵主对角线元素和inv(A)求逆矩阵建议优先改用A\bcond(A)计算条件数判断线性方程组是否病态使用这些函数时最有效的学习方法是在命令窗口输入doc 函数名例如doc norm。文档页面会给出定义、参数说明、算法背景和示例比死记硬背更可靠。注意不要用det(A)是否接近 0 来判断矩阵是否可逆。浮点误差会让行列式结果在极小值和真实值之间产生偏差更可靠的是rcond(A)或cond(A)。3. 解线性方程组用 A\b 而不是盲目求逆3.1 从数学式 Axb 到代码 A\b线性方程组是线性代数最直接的工程应用。数学上求解Ax b需要找到向量x使得矩阵A乘上x得到b。在 MATLAB 中最推荐的写法是A [2 1; 1 3]; b [5; 6]; x A\b;这里的反斜杠符号\是左除专门用于求解Ax b。运行后得到的x就是一个满足条件的解。可以通过残差验证residual A*x - b; norm(residual)如果norm(residual)非常小比如1e-14左右说明解是正确的。因为浮点运算存在误差要求结果为严格的 0 并不现实只要残差在机器精度量级即可。必须区分左除和右除。A\b是左除对应A*x bb/A是右除对应x*A b。两者含义不同不能混用。3.2 为什么不要只依赖 inv(A)有些人解方程会写成x inv(A)*b。这种写法在数学上等价于x A^(-1)*b但在数值计算中通常不是首选。inv(A)会先计算完整逆矩阵再用矩阵乘法和b相乘。这个过程计算量更大而且在矩阵接近奇异或条件数很大时会放大数值误差。A\b则会根据A的特点选择合适算法比如 LU 分解、Cholesky 分解、最小二乘等稳定性更好。下面用 Hilbert 矩阵对比两种方式。Hilbert 矩阵是典型的病态矩阵阶数升高后问题非常明显n 10; A hilb(n); b ones(n, 1); x1 A\b; x2 inv(A)*b; residual1 norm(A*x1 - b); residual2 norm(A*x2 - b); disp([residual1 residual2]);在我本机的一次运行中x1的残差通常比x2小一个数量级以上矩阵阶数越高差异越明显。结论是在 MATLAB 中解Ax b优先写A\b不要用inv(A)*b。只有在确实需要得到逆矩阵本身并且矩阵规模很小、条件数很好时才使用inv(A)。3.3 超定方程、欠定方程与稀疏矩阵实际工程里不是所有方程组都有唯一解。超定方程组指方程数大于未知数个数例如用多个观测点拟合一条曲线。A\b在矩阵不是方阵时不会报错而会返回最小二乘意义下的解即让norm(A*x - b)最小的x。例如A [1 1; 1 2; 1 3]; b [2.1; 2.9; 4.2]; x A\b;欠定方程组指未知数个数大于方程数A\b会返回一个具有最小范数的解。这类问题在控制系统分配、信号重构中很常见。处理大型稀疏矩阵时普通矩阵会消耗大量内存。比如一个只在对角线附近有值的带状矩阵如果按全矩阵存储很多零元素会被白白存下来。可以改用稀疏矩阵n 100; A spdiags([-ones(n,1), 4*ones(n,1), -ones(n,1)], [-1 0 1], n, n); b ones(n, 1); x A\b;spdiags按对角线创建稀疏矩阵。A\b对稀疏矩阵仍然是有效的它会自动选择适合稀疏结构的算法。对于更大的矩阵还可以考虑迭代求解器例如pcgx pcg(A, b, 1e-10, 100);这一行表示用共轭梯度法求解目标相对残差为1e-10最大迭代次数为 100。迭代法只适合满足特定条件的矩阵使用前需要确认矩阵是否对称正定。不要为了展示技术而盲目使用迭代法。4. 矩阵分解和特征值看清矩阵的内部结构4.1 特征值、奇异值与矩阵分解的含义矩阵分解的本质是把一个复杂矩阵拆成几个结构简单的矩阵相乘以便看清它内部的变换特性。特征值分解适用于方阵奇异值分解适用于任意矩阵。用通俗的方式理解特征值描述的是一个矩阵在某个特征向量方向上“拉伸或者压缩了多少倍”。如果特征值为 2说明沿对应特征向量方向的向量被放大 2 倍如果特征值为 0说明该方向上存在零空间。奇异值分解则把特征值的思想推广到非方阵例如把一张m×n的灰度图像分解成左奇异向量、奇异值矩阵和右奇异向量常用于图像压缩和降维。LU 分解会把矩阵拆成下三角矩阵和上三角矩阵适合解线性方程组。QR 分解适合最小二乘问题和特征值算法的底层实现。理解这些分解的适用场景比死记分解公式更有价值。4.2 eig、svd、lu、qr 的最小用法特征值分解用eigA [1 2; 2 1]; [V, D] eig(A); lambda diag(D);D是对角矩阵对角线上的元素就是特征值V的每一列是对应的特征向量。验证时可以检查特征值分解的核心等式err norm(A*V - V*D, fro);err很接近 0说明分解正确。奇异值分解用svdA [1 2 3; 4 5 6]; [U, S, V] svd(A); reconstructed U*S*V; err norm(A - reconstructed, fro);这里U是2×2矩阵S是2×3对角矩阵V是3×3矩阵。svd会根据输入矩阵的形状自动调整输出矩阵大小。LU 分解用luA [4 3; 6 3]; [L, U, P] lu(A); reconstructed P * L * U; err norm(A - reconstructed, fro);注意MATLAB 的lu返回的L、U、P满足P*A L*U所以要重建原矩阵A需要用到P * L * U而不是直接L*U。QR 分解用qrA [1 2; 3 4; 5 6]; [Q, R] qr(A); err norm(A - Q*R, fro);Q是正交矩阵R是上三角矩阵。QR 分解在最小二乘问题中很关键数值稳定性通常优于直接构造法方程。4.3 分解结果如何验证很多学习者在运行eig、svd后只看输出结果从不验证。实际上任何矩阵分解都应该用“残差”验证一下。常用验证命令norm(A*V - V*D, fro) norm(A - U*S*V, fro) norm(A - P*L*U, fro) norm(A - Q*R, fro)fro表示 Frobenius 范数相当于把矩阵当成一个长向量计算二范数适合衡量整体误差。残差如果接近1e-14量级说明分解正确如果很大说明代码里的重建顺序或维度理解有问题。此外评估矩阵是否“稳定”也很重要。条件数cond(A)反映的是线性方程组对误差的敏感程度。条件数接近 1说明问题良态条件数很大说明病态微小扰动都会导致结果剧烈变化。A [1 2; 2 4]; cond(A)这个矩阵是奇异的cond(A)会非常大。实际计算时不要只看某个数字要结合问题的物理背景判断哪些误差可以接受。注意svd得到的S可能是矩形对角矩阵重建和切片时要严格匹配维度。写代码时先查看size(U)、size(S)、size(V)再决定如何相乘。5. 线性代数不是孤立知识拟合、图像与仿真的落地点5.1 用线性代数做最小二乘拟合最小二乘拟合是线性代数最直观的工程应用也是从“会解方程组”到“会用方程组解决数据问题”的关键一步。假设有一组散点数据(x_i, y_i)想用二次多项式y p1 p2*x p3*x^2去拟合。如果把每个数据点代入就能得到p1 p2*x1 p3*x1^2 y1 p1 p2*x2 p3*x2^2 y2 ...写成矩阵形式就是A*p y其中A的第一列是 1第二列是x第三列是x.^2。当数据点个数大于 3 时这是一个超定方程组。MATLAB 代码如下rng(0); x (0:0.1:1); y 1 2*x 3*x.^2 0.1*randn(size(x)); A [ones(size(x)), x, x.^2]; p A\y; yFit A*p; plot(x, y, o, x, yFit, -);A\y返回的就是最小二乘解也就是让拟合曲线与原始数据整体误差最小的系数。要检查拟合质量可以计算残差residual y - yFit; norm(residual)使用rng(0)是为了固定随机数种子让每次运行得到相同数据方便复现。如果不固定种子每次randn产生的噪声不同拟合结果也会有细微差异。5.2 图像旋转和缩放背后的矩阵变换图像在 MATLAB 中本身就是矩阵。灰度图是一个m×n的二维矩阵彩色图是一个m×n×3的三维数组。很多图像处理操作比如旋转、缩放、平移底层都离不开线性变换。二维旋转矩阵是theta pi/6; R [cos(theta), -sin(theta); sin(theta), cos(theta)];对一个坐标点point [1; 0]做旋转point [1; 0]; rotated R * point;rotated就是按逆时针旋转 30 度后的新坐标。实际图像旋转中MATLAB 通常用imrotate等封装函数但理解旋转矩阵可以帮助你判断为什么一个角度很小图像边缘会出现锯齿为什么变换之后需要插值为什么坐标变换顺序会影响最终结果。利用 SVD 做图像近似压缩也是很好的练习思路是保留前k个奇异值丢掉后面较小的奇异值% 示意代码需要先准备灰度矩阵 I [U, S, V] svd(I); k 50; Iapprox U(:,1:k) * S(1:k,1:k) * V(:,1:k);这个例子能直观展示“矩阵分解”的价值只保留少数奇异值就能保留图像的主要结构。5.3 从 Simulink 到光学仿真线性代数是公共语言MATLAB 的工程场景非常分散Simulink 建模、光学工具箱、电扫阵列仿真、GRIB 数据读取、图像处理大作业、随机游走模型看起来互不相关底层都会回到矩阵运算。Simulink 里传递的多维信号通常就是矩阵。矩阵维度不一致时会出现连接错误排查时首先要确认每个信号的size。光学仿真中的光线传输、阵列方向图计算往往通过复矩阵乘法和相位叠加完成。图像处理中的卷积、滤波、PCA全部依赖矩阵操作。随机游走模型是一个适合练习“向量化思维”的例子。不要用大量for循环逐点计算可以把多条轨迹同时组织成矩阵rng(0); nSteps 1000; nWalks 10; steps sign(randn(nSteps, nWalks)); positions cumsum(steps); plot(positions);randn(nSteps, nWalks)生成一个1000×10的矩阵sign把值映射为1或-1cumsum按第一维累积求和相当于同时计算 10 条随机游走轨迹。这种写法简洁也更容易扩展到大规模模拟。6. 常见报错、浮点误差与排错链路6.1 典型报错与排查表MATLAB 线性代数代码的报错大多数集中在维度不匹配、运算符误用、矩阵奇异、函数不存在这几类。实际排错时不要只看报错最后一行要从第一行错误提示开始看因为它往往指向真正出问题的位置。报错或警告常见原因检查方式处理建议Matrix dimensions must agree.矩阵乘法或加减法的维度不匹配size(A)、size(B)对比维度统一内维必要时转置或reshapeDimensions of arrays being concatenated are not consistent.行拼接或列拼接时行列数不一致检查[A B]或[A; B]的尺寸保证横向拼接行数相同纵向拼接列数相同Matrix must be square.对非方阵使用eig、det等函数size(A)查看是否为方阵方阵分解用eig非方阵用svdWarning: Matrix is singular to working precision.矩阵奇异或接近奇异rank(A)、cond(A)、rcond(A)检查数据来源改用pinv或正则化Unrecognized function or variable xxx.函数名拼写错误、工具箱缺失或未定义变量which xxx、ver查看工具箱检查拼写和工具箱安装情况Array indices must be positive integers or logical values.用小数、0 或负值作为索引打印索引变量的具体值检查索引是否由find或计算产生错误6.2 浮点误差、条件数与病态矩阵MATLAB 的数值计算基于浮点数浮点数不能精确表示所有实数。因此0.1 0.2可能不等于精确的0.3矩阵求逆的结果也不是数学意义上的绝对精确值。这不是 MATLAB 的 bug而是所有浮点计算共有的特性。病态矩阵会让浮点误差被放大。典型例子是 Hilbert 矩阵n 8; A hilb(n); cond(A) det(A) rank(A)cond(A)会非常大说明矩阵对误差极其敏感。det(A)虽然很小但行列式小不等于矩阵不可逆。实际判断一个矩阵是否奇异或接近奇异更推荐看条件数和非零奇异值的个数。遇到病态矩阵时常见处理思路包括检查数据是否经过归一化或缩放。优先使用A\b而不是inv(A)*b。尝试pinv(A)求伪逆得到最小范数解。对问题增加正则化项比如岭回归中的(A*A lambda*I)。改用更稳定的分解例如 QR 分解或 SVD。6.3 一套可复用的排错流程线性代数代码出问题时不要直接重新跑一遍建议按下面这条流程排查确认输入数据打印size(A)、size(b)确认矩阵不是空矩阵。确认运算符是矩阵乘法*还是逐元素乘法.*是共轭转置还是普通转置.。确认函数是否存在运行which或doc查看函数来源。确认矩阵状态运行cond(A)、rank(A)判断是否病态或奇异。确认结果是否有异常检查结果中是否出现NaN或Inf并用norm(A*x - b)验证残差。确认算法是否适合判断方程组是超定、欠定还是方阵是否需要稀疏存储或迭代法。可以把这个流程写成一个注释模板放在自己的脚本开头每次排错都按顺序过一遍% 1. size % 2. 运算符 % 3. which % 4. cond / rank % 5. norm residual % 6. 算法选型7. 学习环境与生产实践中的最佳做法7.1 学习环境与生产环境的差异学习阶段主要目标是跑通代码、理解概念允许脚本写得随意。进入项目或生产环境后代码要考虑可维护性、可测试性和性能。两者之间的差异可以用下面这张表说明维度学习环境生产环境目标理解原理、跑通结果稳定运行、可维护、可排错代码组织一段脚本函数 测试脚本输入检查通常省略必须校验维度、类型、取值范围错误处理直接看报错信息捕获异常、记录日志、返回明确提示性能数据规模小怎么写都行向量化、稀疏存储、避免不必要的复制可复现性随机数据即可固定随机种子、保存数据和参数学习时可以只用命令窗口或单个脚本生产环境建议把线性代数计算封装成函数保留输入输出接口。7.2 工程化建议脚本、函数、单元测试一个典型的线性方程组求解函数可以这样组织function x solveLinearSystem(A, b) validateattributes(A, {numeric}, {2d, square}, solveLinearSystem, A, 1); validateattributes(b, {numeric}, {column}, solveLinearSystem, b, 2); if size(A, 1) ~ size(b, 1) error(solveLinearSystem:dimMismatch, A 的行数必须等于 b 的长度); end x A\b; endvalidateattributes负责检查输入类型和维度避免函数内部在运行到一半时才发现错误。列向量检查可以用size(b, 2) 1也可以使用validateattributes的column属性。配合测试脚本A [2 1; 1 3]; b [5; 6]; x solveLinearSystem(A, b); assert(norm(A*x - b) 1e-10, 残差过大);这里assert会在线性方程组的解不满足要求时抛出异常。把测试脚本保存下来后续修改函数时再运行一次能快速发现回归问题。生产环境还需要考虑日志、权限、监控、回滚和异常处理。如果是把 MATLAB 程序部署给其他人使用还要确认目标机器的 MATLAB 版本、工具箱和许可证是否匹配。7.3 后续学习路线建议如果你刚学完这一篇建议按以下顺序继续练习构造一个3×3随机矩阵计算det、rank、cond用A\b求解并验证残差。生成一组带噪声的散点用多项式拟合分别尝试 1 次、2 次、5 次曲线观察过拟合现象。用svd对一张灰度图做低秩近似观察k10、k50、k100时图像的变化。对一个对称矩阵运行eig验证特征向量是否正交。把一个用循环写的随机游走模型改成矩阵版本对比运行时间。理解 MATLAB 线性代数的关键不是记住所有函数而是建立“先看维度、再选算法、最后验证残差”的反射。建议先从A\b开始练跑通后再研究eig和svd。遇到报错时回看本文第 6 节这条路径基本能覆盖 MATLAB 线性代数的大部分日常使用场景。
返回列表