MATLAB向量化编程:利用隐式扩展高效生成矩阵
1. 项目概述告别循环拥抱向量化在MATLAB的日常编程中我们经常遇到一个经典任务如何高效地生成一个矩阵其元素值由行列索引的某种函数关系决定新手甚至一些有经验的开发者第一反应往往是使用双重for循环逐个元素计算并赋值。这很直观但性能上往往不尽如人意尤其是在矩阵维度较大时。循环带来的开销会成为明显的瓶颈。今天要聊的这个技巧核心就是利用MATLAB R2016b版本引入的一项关键特性——算术扩展来彻底摒弃循环以一种更“MATLAB”的方式即向量化操作来优雅地构造矩阵。这不仅仅是写代码风格的变化更是性能上的巨大提升。想象一下你不再需要写for i 1:m, for j 1:n, A(i,j) f(i,j); end, end这样冗长的嵌套而是用一行简洁的表达式直接得到结果矩阵。这行代码不仅运行更快而且更清晰地表达了你的数学意图。这个技巧适合所有使用MATLAB进行科学计算、数据分析或算法开发的工程师和研究人员。无论你是要生成一个坐标网格、一个自定义的滤波器核还是一个依赖于索引的权重矩阵掌握这个方法都能让你的代码更高效、更专业。2. 核心原理算术扩展与隐式扩展要理解这个技巧必须先搞清楚两个核心概念算术扩展和MATLAB R2016b带来的隐式扩展。很多人容易混淆它们但其实它们是紧密相关的两个层面。2.1 算术扩展从标量到矩阵的桥梁算术扩展更广义地说是MATLAB中许多运算符和函数能够自动处理不同尺寸数据的能力。最常见的例子就是标量与矩阵的运算。当你写A 5时MATLAB会自动将标量5扩展为与矩阵A同尺寸的矩阵其所有元素都是5然后执行逐元素加法。这本质上是一种“广播”机制。在R2016b之前这种扩展主要发生在标量与矩阵之间或者两个维度完全相同的矩阵之间。2.2 R2016b的隐式扩展维度的自动对齐R2016b版本的更新将这种“广播”能力推向了一个新的高度即隐式扩展。它允许在两个数组的维度满足特定条件时自动将维度为1的维度进行扩展以匹配另一个数组的对应维度。具体规则是对于两个数组A和B当它们的每个维度大小要么相同要么其中一个为1时这两个数组就是兼容的。MATLAB会自动将维度为1的那个维度复制扩展到与另一个数组对应维度相同的大小然后执行逐元素运算。举个例子列向量v [1; 2; 3](尺寸 3x1)行向量h [4, 5, 6](尺寸 1x3)在R2016b之前v h会报错因为维度不匹配。但在R2016b及之后MATLAB会进行隐式扩展v的第二个维度是1h的第二个维度是3。v的第二维被扩展为3变成[1,1,1; 2,2,2; 3,3,3]逻辑上。h的第一个维度是1v的第一个维度是3。h的第一维被扩展为3变成[4,5,6; 4,5,6; 4,5,6]逻辑上。现在两个操作数都变成了3x3的矩阵可以执行逐元素相加结果是一个3x3矩阵[14, 15, 16; 24, 25, 26; 34, 35, 36] [5, 6, 7; 6, 7, 8; 7, 8, 9]这里有一个至关重要的理解这个扩展过程是在内存中“逻辑上”完成的MATLAB的底层优化确保了它不会真的去创建两个完整的3x3临时矩阵从而在大多数情况下非常高效。2.3 生成矩阵的核心思路我们的目标是用行列索引的函数f(i,j)来填充矩阵A。利用隐式扩展我们可以将行索引i构造成一个列向量将列索引j构造成一个行向量。当这两个向量进行某种运算时隐式扩展就会自动生成一个矩阵其(i,j)位置的元素正好是f(i, j)。假设我们要生成一个m x n的矩阵其中A(i,j) i^2 j。生成行索引向量row_vec (1:m).这是一个m x 1的列向量内容是[1; 2; ...; m]。生成列索引向量col_vec 1:n这是一个1 x n的行向量内容是[1, 2, ..., n]。利用隐式扩展计算A row_vec.^2 col_vec。row_vec.^2得到m x 1的列向量[1; 4; ...; m^2]。col_vec是1 x n的行向量[1, 2, ..., n]。执行加法时row_vec.^2的第二维值为1扩展n次col_vec的第一维值为1扩展m次。结果A(i,j) i^2 j。就这样没有使用一个循环我们就得到了目标矩阵。这个方法的性能远超循环尤其是当m和n很大时。注意在R2016b之前要实现类似功能通常需要借助repmat函数显式地复制向量或者使用bsxfun函数。隐式扩展的引入使得语法更加简洁直观可以看作是bsxfun的语法糖。但性能上三者隐式扩展、bsxfun、优化后的repmat在底层实现上已经非常接近。3. 实战演练从简单到复杂的矩阵构造理解了原理我们通过几个具体例子来看看如何应用这个技巧解决实际问题。我会从最简单的线性网格开始逐步过渡到更复杂的函数关系。3.1 案例一生成二维坐标网格这是图像处理、数值计算中最常见的需求之一。我们需要生成两个矩阵X和Y其中X的每一行都是横坐标向量Y的每一列都是纵坐标向量。假设我们的坐标范围是x从-1到1共Nx个点y从-2到2共Ny个点。传统循环方法低效Nx 100; Ny 150; x linspace(-1, 1, Nx); y linspace(-2, 2, Ny); X zeros(Ny, Nx); Y zeros(Ny, Nx); for i 1:Ny for j 1:Nx X(i, j) x(j); % 注意索引X的列由x决定 Y(i, j) y(i); % Y的行由y决定 end end使用隐式扩展的向量化方法Nx 100; Ny 150; x linspace(-1, 1, Nx); % 1 x Nx 行向量 y linspace(-2, 2, Ny); % Ny x 1 列向量 X x; % 行向量 Y y; % 列向量 % 利用隐式扩展 X repmat(x, Ny, 1); % 方法A显式复制R2016b前常用 % 或者更直接地利用赋值时的扩展但这里需要构造 % 实际上最优雅的方式是 X x; % 先定义为行向量 Y y.; % 转换为列向量 % 然后利用隐式扩展进行运算生成网格时会自动扩展。 % 但为了直接得到X和Y矩阵我们可以这样做 X ones(Ny, 1) * x; % 外积方式但不够“隐式” % 最佳实践利用隐式扩展 X x; % 1 x Nx Y y.; % Ny x 1 % 当我们需要计算每个点的半径 r sqrt(X.^2 Y.^2) 时直接写 % R sqrt(X.^2 Y.^2); % 这里会报错因为X和Y尺寸不直接匹配 % 正确生成网格矩阵的方法是 X linspace(-1, 1, Nx); Y linspace(-2, 2, Ny).; [X_grid, Y_grid] meshgrid(x, y); % meshgrid是专门函数 % 但如果我们想用隐式扩展“模仿”meshgrid X_grid x; % 1 x Nx Y_grid y.; % Ny x 1 % 在计算时例如 Z sin(X_grid Y_grid)MATLAB会自动扩展。 % 但为了得到网格矩阵本身我们需要一个触发扩展的操作 X_matrix ones(size(y.)) * x; % 等价于 repmat(x, length(y), 1) Y_matrix y. * ones(size(x)); % 等价于 repmat(y., 1, length(x))上面的例子展示了思维过程。实际上对于纯网格生成meshgrid或ndgrid函数仍然是首选因为它们语义最清晰。但理解其背后的原理就是隐式扩展或类似机制对于掌握向量化思维至关重要。更地道的“隐式扩展”生成网格方式通常是在需要用到网格的场合直接计算而不是先生成网格矩阵。例如计算函数f(x,y) sin(x) * cos(y)在网格上的值x linspace(0, 2*pi, 100); % 1x100 y linspace(0, pi, 80).; % 80x1 % 直接计算隐式扩展在背后工作 Z sin(x) .* cos(y); % sin(x)是1x100cos(y)是80x1结果Z是80x100这里我们并没有显式生成X和Y矩阵而是在计算表达式中直接使用了行向量x和列向量y。sin(x)和cos(y)的结果仍然保持原来的维度在进行点乘.*时隐式扩展发生直接得到了结果矩阵Z。这才是隐式扩展技巧的精髓让数据保持在合适的维度上让运算规则自动处理扩展。3.2 案例二构造Toeplitz矩阵或类似带状矩阵Toeplitz矩阵是每条对角线元素都相等的矩阵在信号处理中很常见。生成一个Toeplitz矩阵第一列是c第一行是r。传统方法可能会用循环填充对角线。向量化方法可以利用隐式扩展通过索引的差值或和值来生成。c [1, 2, 3, 4].; % 4x1 列向量第一列 r [1, 5, 6, 7]; % 1x4 行向量第一行 (r(1)应与c(1)相同) n length(c); m length(r); % 生成行索引矩阵 I 和列索引矩阵 J % I(i,j) i, J(i,j) j I (1:n).; % n x 1 J 1:m; % 1 x m % Toeplitz矩阵的 (i,j) 元素是 c(i-j1)? 需要仔细定义。 % 更通用的方法是利用差值索引。但标准Toeplitz可以用toeplitz(c,r)函数。 % 这里我们用隐式扩展演示一个类似原理生成矩阵A其中A(i,j) |i - j|。 I (1:n).; % n x 1 列向量 J 1:m; % 1 x m 行向量 A abs(I - J); % 隐式扩展I扩展列J扩展行。 % 结果A就是一个n x m的矩阵A(i,j) |i-j|。这个abs(I - J)就是一个完美的例子。I是n x 1J是1 x m相减时I的每一列都复制m份减去整个J行向量J的每一行都复制n份减去整个I列向量最终得到一个n x m的差值矩阵再取绝对值。3.3 案例三生成自定义核函数矩阵在机器学习或图像处理中我们可能需要一个高斯径向基函数矩阵其中元素K(i,j) exp(-gamma * ||x_i - x_j||^2)x_i和x_j是数据点。假设X是一个m x d的数据矩阵m个样本d维特征。我们需要计算一个m x m的核矩阵K。循环方法极其低效双重循环计算每对样本的距离。向量化方法利用隐式扩展 这个例子稍微复杂需要分解步骤。核心思想是避免计算每对向量差值的显式循环。% 假设我们有数据矩阵 X (m x d) m 1000; d 10; X randn(m, d); gamma 0.1; % 计算核矩阵 K其中 K(i,j) exp(-gamma * ||X(i,:) - X(j,:)||^2) % 技巧利用 ||a-b||^2 ||a||^2 ||b||^2 - 2*a*b % 1. 计算每个样本的范数平方 (m x 1) X_norm2 sum(X.^2, 2); % 沿第二维求和得到 m x 1 列向量 % 2. 计算 Gram 矩阵 (点积矩阵) G X * X (m x m) G X * X; % 3. 利用隐式扩展计算距离平方矩阵 D2 % D2(i,j) X_norm2(i) X_norm2(j) - 2*G(i,j) % X_norm2 是 m x 1 X_norm2 是 1 x m D2 X_norm2 X_norm2 - 2 * G; % 这里发生了关键的隐式扩展 % X_norm2 (m x 1) 加上 X_norm2 (1 x m)扩展为 m x m 矩阵其(i,j)元素为 X_norm2(i) X_norm2(j) % 然后减去 2*G (m x m) % 4. 计算核矩阵 K exp(-gamma * D2);这段代码完全没有使用循环就计算出了所有样本对之间的高斯核函数值。其中X_norm2 X_norm2这一步就是隐式扩展的典型应用一个列向量和一个行向量相加直接生成一个矩阵。这个技巧将计算复杂度从循环的O(m^2 * d)降低到了矩阵运算的O(m^2 * d)但常数项和实际执行效率远优于循环并且代码简洁明了。实操心得在处理这类需要计算所有两两组合的问题时思考如何将计算转化为“一个向量/矩阵与另一个向量/矩阵的运算”是应用隐式扩展的关键。通常这涉及到将问题分解为“自身项”、“交互项”等部分并利用sum(X.^2,2)、X * X等操作先计算好中间量。4. 性能对比与深度解析说一千道一万向量化到底能快多少我们来做一个直观的对比并深入分析其背后的原因。4.1 基准测试循环 vs. 隐式扩展我们以生成一个矩阵A(i,j) sqrt(i) log(j)为例矩阵大小设为2000 x 2000。m 2000; n 2000; % 方法1双重 for 循环 tic; A_loop zeros(m, n); for i 1:m for j 1:n A_loop(i, j) sqrt(i) log(j); end end time_loop toc; fprintf(双重循环耗时: %.4f 秒\n, time_loop); % 方法2使用隐式扩展 tic; I (1:m).; J 1:n; A_vectorized sqrt(I) log(J); % 关键的一行 time_vec toc; fprintf(隐式扩展耗时: %.4f 秒\n, time_vec); % 验证结果是否一致 fprintf(结果最大差异: %e\n, max(abs(A_loop(:) - A_vectorized(:))));在我的测试环境MATLAB R2023a下结果可能是这样的双重循环耗时: 2.8571 秒 隐式扩展耗时: 0.1024 秒 结果最大差异: 0.000000e00性能提升超过20倍这个差距随着矩阵尺寸的增大会更加显著。对于5000 x 5000的矩阵循环方法可能需要几十秒甚至几分钟而向量化方法可能仍然在1秒以内。4.2 为什么向量化如此之快解释器开销MATLAB是解释型语言。每次执行循环体中的语句解释器都需要进行解析、类型检查、函数查找等操作。对于几百万次的循环迭代这个开销是巨大的。而向量化操作将整个计算打包交给底层高度优化的、用C/C或Fortran编写的库函数如BLAS, LAPACK, Intel MKL去执行解释器开销几乎可以忽略不计。内存访问与缓存优化现代CPU具有多级缓存。循环操作往往导致非连续的内存访问模式容易引起缓存失效。而向量化运算底层库函数被设计为以最有利于CPU缓存的方式如循环展开、数据分块来访问内存极大地提高了数据吞吐量。并行计算潜力许多底层的矩阵运算库本身就支持多线程并行。当你执行sqrt(I) log(J)时MATLAB可能会自动调用多线程来并行计算sqrt(I)和log(J)以及最终的加法。这是你“免费”获得的并行加速无需自己写parfor。循环优化JIT的局限MATLAB的即时编译器JIT可以加速简单的循环但对于内部有函数调用如sqrt,log、条件判断或复杂索引的循环JIT的优化效果有限。向量化则完全避开了这些限制。注意事项隐式扩展虽然高效但并非没有成本。它确实会在逻辑上创建扩展后的数组。对于非常大的数组即使底层有优化扩展操作本身也可能消耗可观的内存和计算资源。在极端情况下如果两个输入数组都需要向高维扩展可能会触发巨大的临时内存分配。例如一个10000x1的向量和一个1x10000的向量相加逻辑上会产生一个10000x10000的临时矩阵。虽然MATLAB的优化可能避免完全物化它但在编写代码时仍需对数据规模有意识。如果遇到内存不足的问题可能需要考虑分块计算或使用bsxfun在旧版本中的惰性求值特性。4.3 与repmat和bsxfun的历史对比在R2016b之前要实现类似功能主要有两种方式repmat函数显式复制数组以达到维度匹配。A sqrt(repmat((1:m)., 1, n)) log(repmat(1:n, m, 1));这种方法语法清晰但repmat需要实际创建数据的副本在旧版本中可能带来额外的内存分配和复制开销。bsxfun函数二进制单例扩展函数。这是R2016b前推荐的向量化利器。A bsxfun(plus, sqrt((1:m).), log(1:n)); % 或者对于更复杂的函数句柄 A bsxfun((i,j) sqrt(i)log(j), (1:m)., 1:n); % 注意bsxfun不支持这种匿名函数直接用于两个数组这里仅为示意。实际需分开计算。 % 正确用法是 I (1:m).; J 1:n; A bsxfun(plus, sqrt(I), log(J));bsxfun在底层实现了惰性扩展避免了repmat可能带来的额外内存占用性能通常更好。R2016b隐式扩展的引入本质上将bsxfun的功能集成到了基本的算术运算符,-,.*,./,.^,,等和一部分函数中。现在A B在维度兼容时其行为就等同于bsxfun(plus, A, B)。这使得代码更加简洁易读。兼容性考虑如果你的代码需要在R2016b之前的版本中运行那么必须使用bsxfun或repmat。对于新项目或维护中的项目如果确定环境在R2016b以上应优先使用隐式扩展因为它语法最自然。5. 高级技巧与边界情况处理掌握了基础用法后我们来看看一些更高级的场景和需要注意的坑。5.1 处理多维数组N维数组隐式扩展不仅限于二维矩阵对N维数组同样有效。规则是一样的对应维度大小要么相等要么其中一个为1。例如我们有一个3维数组A尺寸a x b x c和一个2维数组B尺寸1 x b x c那么A B是合法的B会在第一个维度上扩展a次。假设我们要对一个3D体数据例如128x128x64的MRI图像的每一帧128x128加上一个不同的2D背景场1x128x128。volume randn(128, 128, 64); % 3D 体数据 background_field rand(1, 128, 128); % 2D背景场但定义为 1x128x128 % 我们想将 background_field 加到 volume 的每一帧上 corrected_volume volume background_field; % 隐式扩展发生background_field 的第一维1扩展为128这里background_field的维度是[1, 128, 128]volume的维度是[128, 128, 64]。检查对应维度维1: 1 vs 128 - 兼容1扩展为128维2: 128 vs 128 - 兼容相等维3: 128 vs 64 -不兼容两者都不为1且不相等。 所以这个操作会报错“数组维度必须匹配或者为单例维度可扩展”。正确的做法是背景场应该与每一帧的xy平面对应其维度应为[128, 128, 1]。background_field_2d rand(128, 128); % 2D矩阵128x128 % 为了与volume相加我们需要将其变为 128x128x1 background_field_3d background_field_2d; % MATLAB中2D数组可视为第三维为1 % 或者显式 reshape % background_field_3d reshape(background_field_2d, [128, 128, 1]); corrected_volume volume background_field_3d; % 现在兼容了128x128x1 扩展第三维到64这个例子提醒我们在使用隐式扩展时必须非常清楚每个数组的维度并确保它们按照规则是兼容的。permute、reshape、squeeze等函数是调整维度的好帮手。5.2 与逻辑索引和关系运算的结合隐式扩展同样适用于关系运算符,,,~等和逻辑运算符,|,~。这可以用于生成复杂的逻辑掩膜。例如我们想从一个矩阵中筛选出那些值大于其所在行平均值同时小于其所在列平均值的元素。A rand(5, 7); row_mean mean(A, 2); % 5x1每行的均值 col_mean mean(A, 1); % 1x7每列的均值 % 生成逻辑掩膜 mask (A row_mean) (A col_mean); % 隐式扩展在 和 比较中发生 % A row_mean: 比较 5x7 矩阵和 5x1 列向量列向量第二维扩展为7 % A col_mean: 比较 5x7 矩阵和 1x7 行向量行向量第一维扩展为5 % 最终 mask 是一个 5x7 的逻辑矩阵 result A(mask); % 获取满足条件的元素一行代码就完成了复杂的条件判断逻辑非常清晰。5.3 性能陷阱无意中的巨大临时数组虽然隐式扩展很高效但如果不小心也可能写出低效甚至内存爆炸的代码。关键在于意识到哪些操作会触发扩展。反面例子% 假设我们有一个很大的列向量 data (n x 1)和一个权重行向量 weights (1 x m) data randn(1000000, 1); % 100万行 weights linspace(0, 1, 1000); % 1000个权重 % 目标计算加权后的矩阵每列是 data 乘以一个权重 % 错误写法可能内存爆炸 weighted_matrix data * weights; % 这是矩阵乘法结果是 1000000 x 1000 的矩阵占内存约 8GB。 % 这其实不是隐式扩展是矩阵乘法。但我们的本意可能是逐元素缩放。 % 如果本意是每列乘以不同的权重因子应该用 .* % 但 data 和 weights 维度不直接匹配。如果我们的本意是生成一个矩阵M其中M(i,j) data(i) * weights(j)那么data * weights这个外积计算是正确的但它会显式生成一个巨大的稠密矩阵100万 x 1000双精度约8GB这可能超出内存。更优化的思路如果后续操作是求和或其它可以分解的运算也许可以避免生成完整矩阵。例如计算加权和% 计算 sum_{i,j} data(i) * weights(j) * some_kernel(i, j) % 如果 some_kernel 可分离也许能优化。 % 如果必须生成矩阵且内存不足则需要分块处理。**正确且高效的写法用于生成矩阵**就是利用隐式扩展weighted_matrix data .* weights; % 这会报错维度不匹配 % 正确写法 weighted_matrix data .* weights; % 仍然不对 % 需要确保维度兼容。我们需要 data 是 n x 1, weights 是 1 x m weighted_matrix data .* weights; % data (n x 1) .* weights (1 x m) 是合法的隐式扩展 % 但注意.*是逐元素乘法要求维度完全匹配或可扩展。 % 这里 data 是 n x 1, weights 是 1 x m扩展后得到 n x m 矩阵正是我们想要的。 % 所以正确的代码是 weighted_matrix data .* weights; % 前提是 weights 是行向量 % 或者更安全地 weighted_matrix data .* reshape(weights, 1, []); % 确保 weights 是行向量这个操作在逻辑上会产生一个巨大的矩阵但MATLAB的执行是高效的。然而如果n和m都非常大你仍然需要警惕内存消耗。在这种情况下可能需要重新评估算法是否真的需要物化这个完整的矩阵或者是否可以流式处理/分块处理。5.4 调试技巧如何检查扩展是否按预期工作当你对隐式扩展的结果不确定时可以用小数据测试。使用size函数在操作前后检查数组的维度。构造微型测试用例用2x1和1x3这样的向量测试手动计算验证。利用repmat模拟如果不确定隐式扩展的行为可以先用repmat显式扩展看看结果是否一致。这有助于理解扩展逻辑。A [1; 2]; % 2x1 B [3, 4, 5]; % 1x3 C_implicit A B; % 模拟 A_rep repmat(A, 1, size(B,2)); % 将A的列复制3次 - 2x3 B_rep repmat(B, size(A,1), 1); % 将B的行复制2次 - 2x3 C_explicit A_rep B_rep; isequal(C_implicit, C_explicit) % 应该返回 true6. 常见问题与解决方案速查在实际使用中你可能会遇到以下问题。这里提供一个快速排查指南。问题现象可能原因解决方案错误数组维度必须匹配或者为单例维度可扩展。参与运算的两个数组在某个维度上大小既不相等也不为1。1. 使用size(A)和size(B)检查两个数组的维度。2. 使用permute、reshape或squeeze调整其中一个数组的维度使对应维度满足兼容规则。3. 如果其中一个数组是行向量或列向量确保其方向正确。例如需要列向量时用(:)或.转换。代码在R2016a或更早版本中报错旧版本MATLAB不支持算术运算符的隐式扩展。1. 将A B之类的操作替换为bsxfun(plus, A, B)。2. 或者使用repmat显式扩展数组。内存使用量激增甚至“内存不足”隐式扩展逻辑上创建了非常大的中间数组。例如将10000x1向量与1x10000向量相加得到10000x10000矩阵。1.重新思考算法是否真的需要生成这个完整的大矩阵能否用循环分块处理结果是否可以直接计算而不物化中间矩阵2.使用bsxfun在旧版本中bsxfun有时能更高效地处理内存但R2016b后隐式扩展已优化。3.使用稀疏矩阵如果结果矩阵是稀疏的考虑使用稀疏格式存储。4.增加物理内存或使用磁盘交换最后的手段。结果与循环计算结果有微小差异浮点数计算顺序不同导致的舍入误差。循环是逐元素顺序计算向量化是批量计算计算顺序可能不同。这是浮点数运算的固有特性通常差异在eps机器精度量级。如果应用对精度极其敏感需要分析算法稳定性或考虑使用高精度计算工具箱。对于绝大多数工程和科学应用这种差异可以忽略。逻辑运算结果不符合预期可能混淆了元素级逻辑运算, 和短路逻辑运算,想要对更高维数组的特定维度进行操作隐式扩展规则是固定的有时需要手动调整维度。结合permute函数交换维度使你想要广播的维度变为1。例如对3D数组A (a x b x c)想对每个a-b平面加上一个值向量v (b x 1)可以v_perm permute(v, [3,1,2]); % 变成 1 x b x 1result A v_perm; % v_perm在第一维和第三维扩展最后再分享一个我个人的小技巧当你构思一个向量化操作时先在纸上或脑子里画出输入数组的“形状”维度然后像玩拼图一样思考如何通过转置、重塑或增加单例维度让它们的形状变得“兼容”。这常常能帮你快速找到正确的向量化路径。例如面对“对矩阵的每一行减去该行的均值”这种问题立刻想到行均值是列向量(m x 1)原矩阵是(m x n)列向量的第二维是1可以扩展所以直接A - mean(A,2)即可。多练习这种思维会变成你的第二本能。