MATLAB小贴士(1)

这学期担任《系统工程导论》课的助教，看到同学们的作业，看到大家作业里各种青涩可爱的代码，看到大家说自己的 MATLAB 都是从自控课和系统工程课上学的，不由得想为大家做一点事，希望大家在学习使用 MATLAB 的过程中少走一点弯路，希望能对大家以后的学习和科研带来一点帮助。 MAT 的 LAB 我们从进入大学一开始学习了 C 语言，后来又学习了 C++，也许有同学还接触到了 Java，但是 MATLAB 和他们都有点不一样。从名字就可以看出来了，MATLAB 是针对矩阵运算的，所以在 MATLAB 里最欢迎的就是矩阵。MathWorks 里绝顶聪明的工程师们利用各种令人叹服的手段，榨干了计算机的运算能力，将 MATLAB 中的矩阵运算做到世界第一。不相信的同学可以试试看，自己写一段 C 代码实现矩阵乘法，和 MATLAB 的矩阵乘法相比较，看看谁更快。 既然 MATLAB 针对矩阵运算做了大量的优化，那么我们在用 MATLAB 的时候就尽量利用这种优势，将我们的问题用矩阵进行表示，再用矩阵运算来解决。 什么意思呢？ 举个例子吧。比如有一个 n*n 的矩阵 A，我想把第1列乘以1，第2列乘以2，第3列乘以3...怎么实现呢？比如用 C 语言的风格可能会写成这样： for i = 1 : n
   for j = 1 : n
      A(j,i) = A(j,i) * i;
   end
end 在我的 MATLAB 上，对一个 5000*5000 的矩阵，以上代码运行时间为 0.65s 那有没有一种「更矩阵」的方式来达到同样的目的呢？当然可以，我们知道列向量和一个标量可以进行数乘，所以可以这样写： for i = 1 : n
   A(:,i) = A(:,i) * i;
end 上面A(:,i)表示矩阵 A 的第 i 列，类似的A(i,:)表示矩阵 A 的第 i 行。在我的 MATLAB 上，对同样的矩阵只要 0.20s 那能不能更快呢？如果有一个对角矩阵 D 他的对角线元素分别是从 1 到 n，那么这个矩阵和 A 相乘，就可以达到我们的目的。于是这段程序还可以这样写： A = A * sparse(1:n, 1:n, 1:n); 一句话就搞定了，这里sparse(1:n, 1:n, 1:n)这条语句产生了一个对角线元素从 1 到 n 的稀疏矩阵。在我的 MATLAB 上，这句话只花了 0.14s 好了，这个例子很简单，但从中也可以看出，不同的编码方式，对程序运行效率的影响是很大的。最后这种方法和第一种方法相比，效率提高了接近 5 倍。同学们在大一刚入学的时候就学习了 C 语言，所以很多同学仍然将 MATLAB 当做另一个 C 来使用，这并不是很好的使用 MATLAB 的方式。在 MATLAB 里尽量避免for循环，尽量利用矩阵的运算，或者 MATLAB 内置的一些函数进行计算。接下来我将大家在作业中表现出来的问题简单地分析一下，看看都可以做出哪些改进。 对数据进行标准化 有 n 个 m 维的数据点，减去均值后再除以标准差，就得到了均值为 0 方差为 1 的标准化数据。（扯开去说一点，均值的英文是 mean，方差是 variance 常缩写为 var，标准差是 standard deviation 常缩写为 std，起变量名的时候可以参考）这个过程简单吗？当然简单，于是我看到有类似下面这种 C 风格的代码（我整理的，并非出自哪一位同学的作业）： % 计算均值
for i = 1:n
   for j = 1:m
      x_mean(j) = x_mean(j) + x(i,j);
   end
end
x_mean = x_mean / n;

% 计算方差
for i = 1:n
   for j = 1:m
      x_var(j) = x_var(j) + (x(i,j) - x_mean(j))^2;
   end
end
x_var = x_var / (n-1);

% 下面进行标准化
for i = 1:n
   for j = 1:m
      x(i,j) = (x(i,j) - x_mean(j)) / sqrt(x_var(j));
   end
end 这样的代码很清晰易懂，几乎是标准的 C 语言写法。在我们的作业中 n = 3114，m = 14，上面这样的代码在我的 MATLAB 里运行 1000 次需要 5.98s。那么有没有改进的余地呢？前面已经说了，在 MATLAB 里尽量避免使用for循环，尽量使用矩阵运算以及 MATLAB 内置的函数。MATLAB 里有个mean()函数可以计算均值，有个var()函数可以用来计算方差，有个std()函数可以用来计算标准差。好，那么我们就用上这些函数： % 计算均值
x_mean = mean(x);

% 计算标准差
x_std = std(x);

% 下面进行标准化
for i = 1:m
   x(:,i) = (x(:,i) - x_mean(i)) / x_std(i);
end 这段代码运行 1000 次只需要 0.74s，相比原先的代码，效率提高了 8 倍，而且代码更简洁了。 特征值和特征向量 在作业中需要求某个矩阵的特征值和特征向量，当然同学们都会使用 MATLAB 内置的函数[V, D] = eig(A)，这样得到的对角矩阵 D 其对角线元素就是从小到大排列的的特征值。有的同学希望能将这个对角线元素单独抽取出来变成一个行向量或者列向量，还希望特征值从大到小排列，为了达到这个目的，大家用了各种各样的方法。其实 MATLAB 内置的函数就有diag()可以将一个对角矩阵转化为一个向量，也可以将一个向量转化为一个对角矩阵。转换完之后对向量进行逆序排列，我看到有同学硬是写了一个冒泡排序来完成这个工作，你的数据结构老师一定非常欣慰。但是在 MATLAB 里不用从头开始，直接调用sort()函数就行了~但是在这里还有更简单的方法，wrev()函数就可以将一个向量逆序排列，fliplr()函数可以将一个矩阵左右镜像对称，所以如果想得到从大到小排列的特征值，只要这样就行了： % 得到从大到小排列的特征值
[V, D] = eig(A);
lambda = wrev(diag(D));
V = fliplr(V); 这样特征值和特征向量就都排列好了。 在作业中还需要计算相对误差，要将特征值向量从大到小依次累加，直到大于某个阈值。所以可以直接这样写： % 根据相对误差决定主成分个数
lambda_sum = sum(lambda);
tem_sum = 0;
for i = 1:n
   tem_sum = tem_sum + lambda(i);
   if tem_sum / lambda_sum > 1 - rerr
      break;
   end
end 这样当循环结束的时候，i 的值就是需要求的了。但是在 MATLAB 里可以写得更简洁一些： % 根据相对误差决定主成分个数
m = find(cumsum(lambda) / sum(lambda) > 1 - rerr, 1, 'first'); 大家可以利用 MATLAB 的帮助文档，自己查一查上面这句话里各个函数的作用。 另外，如果大家线性代数学得好，直到矩阵的奇异值（Singular Value）的概念，那么大多数情况下使用奇异值分解比使用特征值分解的数值稳定性更好。这一点需要的知识点太多，就不展开说明了。感兴趣的同学可以看看奇异值分解的概念，利用 MATLAB 的帮助查找一下 MATLAB 中奇异值分解的函数是怎么样的。