MATLAB小贴士(2)

没想到上一篇日志有这么多转发...非常感谢各位朋友们的捧场支持~~【鞠躬】这里我还是说明一下，我写这些日志呢，出发点是在批改同学们的作业过程中觉得大家没有用好手中的这个工具，希望能通过一些小贴士，让同学们能够更好地理解 MATLAB 编程的特点，希望能让同学们少走一些弯路。所以我的题目叫「小贴士」而非「教程」，里面的例子也都是同学们作业中遇到的一个一个的小问题，写起来难免有不成体系之感，还请大家原谅。各位朋友们的评论也让我受益匪浅，也欢迎各种讨论和批评~感谢大家！MATLAB 的哲学原谅我很自大地用了个这么闪亮的标题。正如在上一篇日志中提到的，很多初学 MATLAB 的同学，仅仅是把它当做一个简易的、不用定义变量的 C 语言来用，写的代码仍然沿用了 C 代码的思路和风格。这并不是不行，但是（套用一下这个标题）这不是 MATLAB 的哲学，或者说，这不是「工程的思想」（现在洲哥还在用这个梗吗？）那么，MATLAB 的哲学是怎样的呢？用一句话概括就是：相比于具体实现，你更应该关系整体的思路和算法；你所想到的功能，MATLAB 都内置有相应的实现。MATLAB 就像胶水，把各个「小工具」「小函数」黏在一起，造出你想要的轮子。在上一篇日志中提到过，求均值有mean()，求方差有var()，求和有sum()，不止如此，MATLAB 还有各种各样丰富的工具箱，只需用几行语句、几个模块，就能搭建起一个复杂的系统。甚至还能调用 Java、C 等其他语言写的模块。继续看一些具体的例子吧点与点的距离这次作业是 k-means 算法，所以涉及到大量的计算点与点的距离。假如有两个向量集合X和Y，每一行代表一个点，X有 n 行，Y有 m 行，要计算他们之间两两点对的距离，存到一个 n×m 的矩阵里，怎么办呢？这并不困难，对于大部分同学来说很容易想到这样编写：[n, dim] = size(X);  % size() 函数返回一个矩阵的尺寸
m = size(Y, 1);      % size(Y, 1) 表示 Y 的行数
d = zeros(n, m);
for i = 1:n
   for j = 1:m
      for k = 1:dim
         d(i, j) = d(i, j) + (X(i, k) - Y(j, k))^2;
      end
      d(i, j) = sqrt(d(i, j));
   end
end拜 MathWorks 优秀的工程师们所赐，MATLAB 对 for 循环的优化从未间断，在我的版本（2010a）上这段代码说实话跑得并不慢。特别是，如果向量维数不大的话，这段代码还是挺快的。对于 n=3000, m=3 的情况（以下都用这个情况作测试），这段代码只用了 4.7 毫秒。在上一篇日志里我说尽量少用 for 循环，尽量用向量或者矩阵的形式进行计算，在 MATLAB 编程里专门有个词 vectorize 用以描述将数据「向量化」地计算。那么这里我们就可以不针对一个一个的矩阵元素计算，而针对一列一列，或者一行一行进行整体计算，比如可以这样写：for j = 1:m
   for k = 1:dim
      d(:, j) = d(:, j) + (X(:, k) - Y(j, k)).^2;  % 「.^」表示对矩阵每一个元素进行乘方运算
   end
end
d = sqrt(d);  % 最后整体对矩阵 d 开方，避免对 sqrt() 函数的反复调用这段代码对同样的数据只需要 0.35 毫秒，速度提高了 10 倍多。但是！！根据前面提到的 MATLAB 哲学，作为牛逼闪闪的 MATLAB，对于「求两个点集合之间两两元素的距离」这种工程上很常用的功能，怎么可能不提供内置的实现呢？于是我们找到了pdist2(X, Y)这个函数，其作用就是计算X和Y两个集合之间两两点的距离，而且不仅可以算欧氏距离，还可以用其他各种距离，马氏距离、米科夫斯基距离、棋盘格距离……所以我们的代码可以简化成这个样子：d = pdist2(X, Y);一行就搞定了。这行代码在同样的数据上运行需要 0.33 毫秒。在速度上和前一个版本差不多，不过代码的简洁性是无法比拟的。数据归类在 k-means 算法中，计算出距离矩阵后，要将点归入距离中心最近的那个集合，用一个变量label记录这个点的归属。这也是一个很直观的任务，我们可以不假思索地写出代码：for i = 1:n
   label(i) = 1;
   d_min = d(i, 1);
   for j = 1:m
      if d(i, j) < d_min
         d_min = d(i, j);
         label(i) = j;
      end
   end
end对于我们的测试数据，用时 1.15 毫秒。有没有更简洁优雅的方式呢？分析一下这个过程我们可以知道，这无非就是对距离矩阵的每一行进行比较大小，找出最小的那一个的下标。求最小值（顺带找出这个最小值的下标），这也是一个工程上很常用的功能，翻看 MATLAB 的帮助文档可以看到min()函数有一个变种用法是[C, I] = min(A, [], dim)就符合我们的要求：[~, label] = min(d, [], 2);  % 波浪线是占位符，表示忽略这个输出值一行就够了。这段代码运行同样的测试数据只需要 0.067 毫秒，速度提高了差不多 20 倍。在大家提交作业中有一位同学就是这样写的，非常好~他的代码是最短的，同时运行时间也远小于其他同学，比最慢的要快了3000多倍，很不错。数据的选取将数据点分好类之后，需要计算每一类的重心，也就是需要将归属于每一类的点分别选取出来。最直观的代码可以这样写：centroid = zeros(m, 2);     % 记录类别的中心， m 是类别数
cluster_num = zeros(m, 1);  % 记录类别中的点的数目
for i = 1:n
   for j = 1:m
      if label(i) == j
         centroid(j, :) = centroid(j, :) + X(i, :);
         cluster_num(j) = cluster_num(j) + 1;
      end
   end
end
for j = 1:m
   centroid(j, :) = centroid(j, :) / cluster_num(j);
end这段代码耗时 2.66 毫秒。简单分析可以知道，这段代码还是做了一些无用功的，比如最里面的那个判断就无必要，因为label(i)直接指示了对应的类别，因此可以直接根据label(i)的值对相应类别进行操作。可以改进成这样：centroid = zeros(m, 2);     % 记录类别的中心， m 是类别数
cluster_num = zeros(m, 1);  % 记录类别中的点的数目
for i = 1:n
   centroid(label(i), :) = centroid(label(i), :) + X(i, :);
   cluster_num(label(i)) = cluster_num(label(i)) + 1;
end
for j = 1:m
   centroid(j, :) = centroid(j, :) / cluster_num(j);
end减少了很多不必要的判断，这段代码运行时间 1.74 毫秒；能不能更简洁一些呢？我们知道mean()函数可以计算重心（就是均值），但是麻烦在于如何把每一类的数据取出来单独计算。我们知道X(1, :)表示选出了X的第一行，类似的X([1, 3, 5], :)表示选出了X的第 1, 3, 5 行，X(1:100, :)表示选出了X的前 100 行。那么我想选出label变量中等于 1 对应的X中的那些行，怎么办？可以这样写：X(label == 1, :)前提是label的长度必须要和X的行数相等。label == 1这就产生了一个和label一样长的 0/1 向量，放到X的表达式中就表示选出了 1 对应的那些行。所以上面的代码可以写成这样：centroid = zeros(m, 2);     % 记录类别的中心， m 是类别数
for j = 1:m
   centroid(j, :) = mean(X(label == j, :));
end只需要 0.39 毫秒就跑完了。如何找到需要的功能说了这么多，大家也许要问了：作为初学者，我怎么知道哪些功能 MATLAB 有内置的函数，哪些功能要自己手动写呢？很简单，第一步，请坚定地相信上面的 MATLAB 哲学；第二步，打开 Google 输入关键字「MATLAB + 你想要的功能描述」比如上面提到的计算X和Y集合中点与点之间的距离，在 Google 中输入「MATLAB distance between two」这个时候 Google 就自动提示「matlab distance between two sets of points」——可见这是很多人都会遇到的问题，于是确定，看到结果第一条给出的就是 MathWorks 的 MATLAB 在线帮助文档pdist2()函数（当然需要看一下版本号，不同版本的 MATLAB 函数不太一样）使用得多了，自然就知道哪些功能大概是什么函数了。另外，不到万不得已请不要用中文搜，中文资料充斥着大量重复的抄袭和低质量的转载，不一定能快速找到你想要的东西。顺便一说，对于 k-means 这种超级平民的聚类方法，MATLAB 当然也有内置的函数IDX = kmeans(X, k)，这次作业就是让大家自己编写一下这个函数，加深对所学知识的理解。今后在工程上要用到的话，请毫不犹豫地调用 MATLAB 内置的函数吧~再顺便一说，这次作业中程序最快的同学代码的运行时间已经和 MATLAB 内置的kmeans函数运行时间一样了，赞一个！