为什么在 Mathematica 中使用循环是低效的

发布网友 发布时间：2022-04-20 07:06

我来回答

共2个回答

懂视网 时间：2022-05-11 06:59

比如求一个平面稳态导热问题，控制方程就是拉普拉斯方程：

(我才发现原来有[插入公式]这个功能)

按照最简单的毅种循环来写就是：

def laplace(u):
 nx, ny = u.shape
 for i in xrange(1,nx-1):
  for j in xrange(1, ny-1):
  u[i,j] = ((u[i+1, j] + u[i-1, j]) * dy2 + (u[i, j+1] + u[i, j-1]) * dx2) / (2*(dx2+dy2))

你们都不知道numexpr的么←_←
比numpy还黑的科技→_→

虽然能用的运算没多少吧但是对大矩阵的整体运算还是很快的←_← 最近正好在学numpy这个模块。题主可以看看这个教程，不是很全，但是科学计算方面还是有不少东西的：NumPy-快速处理数据
引用教程中的代码：

import time
import math
import numpy as np
x = [i * 0.001 for i in xrange(1000000)] # 初始化数组0.000~999.999 
start = time.clock()
for i, t in enumerate(x):  # 用循环计算正弦值
 x[i] = math.sin(t)  
print "math.sin:", time.clock() - start

x = [i * 0.001 for i in xrange(1000000)]
x = np.array(x)    # 初始化矩阵（这里是一维）
start = time.clock()
np.sin(x,x)    # numpy的广播计算（代替循环）
print "numpy.sin:", time.clock() - start
# 

用numpy, Cython, 或者 weave
Speed up Python
SciPy官网有关于如何提高Python Performance的教程
PerformancePython
用Pyrex/Cython或者weave基本上可以达到C++的速度。
Laplace的例子，500*500矩阵，100次循环。
numpy和pandas.DataFrame的矩阵运算可以广播，可以map。 第一个技巧是，用map和lambda表达式来生成你要的迭代参数，比如生成一个平方表：map(lambda x: x*x, xrange(100))，这是个黑科技，可以很快速的生成你需要的循环参数；
第二个技巧是，熟练使用矩阵掩膜（mask）来简化循环，比如把矩阵a中小于100的值都置零：a[a<100] = 0，比循环快很多；
第三个技巧是，多使用各种库，如numpy, scipy(signal库简直好顶赞)，如果你做图像，opencv库是唯一的选择。
大致是这样，实际应用中更多的是前两个trick混合使用。 想要快，就内嵌C，Python是解释性语言，会比较慢。
有成熟的计算软件时用的C/C+++python的模式，核心算法和耗时最多的逻辑用C/C++,其他用python.

热心网友 时间：2022-05-11 04:07

因为在Mathematica中使用循环确实是低效的。。。。。。
深层次的原因涉及到Mathematica的底层实现所以我不太懂，但是至少从下面几个例子可以看出Mathematica里确实有很多比循环更好的方法

求和
首先举一个最简单的求和例子，求的值。为了测试运行时间取n=10^6
一个刚接触Mathematica的同学多半会这样写

sum = 0;
For[i = 1, i <= 10^6, i++,
sum += Sin[N@i]];
(*其中N@i的作用是把整数i转化为浮点数，类似于C里的double*)
sum

为了便于计时用Mole封装一下，运行时间是2.13秒，如下图

然后一个有一定Mathematica经验的同学多半会知道同样作为循环的Do速度比For快，于是他可能会这么写
然后一个有一定Mathematica经验的同学多半会知道同样作为循环的Do速度比For快，于是他可能会这么写

sum = 0;
Do[sum += Sin[N@i], {i, 1, 10^6}];
sum

如下图，用时1.37秒，比For快了不少

当然了知道Do速度比For快的同学不太可能不知道Sum函数，所以上面其实是我口胡的，他应该会这么写

Sum[Sin[N@i], {i, 1, 10^6}]

如下图，同样的结果，只用了不到0.06秒

如果这位同学还知道Listable属性并且电脑内存不算太小的话，他也可能会这么写

Tr@Sin[N@Range[10^6]]

如下图，只用了不到0.02秒，速度超过For循环的100倍

当然了这只是一个最简单的例子，而且如果数据量更大的话最后一种方法就不能用了。但是这也足以说明在求和时用循环是低效的，无论是内置的Sum函数还是向量化运算，在效率上都远远高于循环
（这部分模仿了不同程序员如何编写阶乘函数这篇文章，强烈推荐对Mathematica有兴趣的同学去看看）

迭代
接下来举一个迭代的例子，（即Logistic map），取，为了测试运行时间同样取n=10^6
还是先用For循环的做法

x = 0.5;
For[i = 1, i <= 10^6, i++,
x = 3.5 x (1 - x);
];
x

如下图，运行时间2.06秒

（Do循环和For类似，篇幅所限这里就不写了，有兴趣的同学可以自行尝试）
（Do循环和For类似，篇幅所限这里就不写了，有兴趣的同学可以自行尝试）
然后看看内置的Nest函数

Nest[3.5 # (1 - #) &, 0.5, 10^6]

如下图，用时0.02秒，又是将近两个数量级的效率差异

当然了Nest的使用涉及到纯函数，对于Mathematica初学者来说可能有些难以理解，而且一些比较复杂的迭代不太容易写成Nest的形式，但是在迭代时Nest（还包括Fold）的效率确实要好于循环
当然了Nest的使用涉及到纯函数，对于Mathematica初学者来说可能有些难以理解，而且一些比较复杂的迭代不太容易写成Nest的形式，但是在迭代时Nest（还包括Fold）的效率确实要好于循环
遍历列表
依然举一个简单的例子：求一个列表中偶数的个数。为测试生成10^6个1到10之间的随机整数

list = RandomInteger[{1, 10}, 10^6];
(*生成10^6个随机整数*)

如果用For循环的话代码是这样的

num = 0；
For[i = 1, i <= 10^6, i++,
If[EvenQ@list[[i]], num++]
];
num

如下图，用时1.73秒

保留上面的思路，单纯的将For循环改为Scan （相当于没有返回结果的Map），代码如下

num = 0；
Scan[If[EvenQ@#, num++] &, list];
num

如下图，用时0.91 秒

（Do循环用时1.00秒左右，篇幅所限就不传图了）

摒弃循环的思路，用其他内置函数写

Count[list, _?EvenQ] // AbsoluteTiming
(*直接用Count数出list中偶数的个数*)
Count[EvenQ /@ list, True] // AbsoluteTiming
(*用Map对list中的每个数判断是否偶数，然后用Count数出结果中True的个数*)
Select[list, EvenQ] // Length // AbsoluteTiming
(*选取list中的所有偶数，然后求结果列表长度*)
Count[EvenQ@list, True] // AbsoluteTiming
(*利用EvenQ的Listable属性直接判断list的每个数是否偶数，然后数出结果中True的个数*)
Sum[Boole[EvenQ@i], {i, list}]
(*对list中的每个元素判断是否偶数，将结果相加*)

结果如下图

这个遍历的例子举得不算特别恰当，但也能说明一些问题了：Mathematica中内置了许多神奇的函数，其中大部分只要使用得当效率都比循环高（而且不是一点半点）。就算非要用循环，也要记得（任何能用Do代替For的时候）
这个遍历的例子举得不算特别恰当，但也能说明一些问题了：Mathematica中内置了许多神奇的函数，其中大部分只要使用得当效率都比循环高（而且不是一点半点）。就算非要用循环，也要记得（任何能用Do代替For的时候）
Do比For快
，（遍历列表时）
Scan比Do快
用向量(矩阵)运算代替循环
这个例子来自如何用 Python 科学计算中的矩阵替代循环？ - Kaiser 的回答，我只是把代码从Python翻译成了Mathematica而已。选这个例子是因为它有比较明确的物理意义，而且效率对比非常明显
代码如下

AbsoluteTiming[
n = 100;
u = unew = SparseArray[{{1, _} -> 1}, {n, n}] // N // Normal;
For[k = 1, k <= 3000, k++,
For[i = 2, i < n, i++,
For[j = 2, j < n, j++,
unew[[i, j]] =
0.25 (u[[i + 1, j]] + u[[i - 1, j]] + u[[i, j + 1]] +
u[[i, j - 1]])
]
];
u = unew;
];
u1 = u;
]
(*用三重循环，迭代3000次*)
ArrayPlot[u1, DataReversed -> True, ColorFunction -> "TemperatureMap"]
(*用ArrayPlot绘图*)

AbsoluteTiming[
n = 100;
u = SparseArray[{{1, _} -> 1}, {n, n}] // N // Normal;
Do[
u[[2 ;; -2, 2 ;; -2]] =
0.25 (u[[3 ;; -1, 2 ;; -2]] + u[[1 ;; -3, 2 ;; -2]] +
u[[2 ;; -2, 3 ;; -1]] + u[[2 ;; -2, 1 ;; -3]]),
{k, 1, 3000}];
u2 = u;
]
(*用矩阵运算，迭代3000次*)
ArrayPlot[u2, DataReversed -> True, ColorFunction -> "TemperatureMap"]
(*用ArrayPlot绘图*)

运行结果For循环用时136秒，矩阵运算用时不足0.5秒，且两者答案完全一样。在算法完全相同的情况下两种写法有着超过200倍的效率差距
（图片太长了这里就不直接显示了，链接放在下面）
http://pic4.mg.com/65d66161f4c674b1149651c32f69935f_b.png

===========================我是结尾的分隔线===============================
这个答案其实从一开始就跑题了，还写了这么长的目的就在于希望让大家切实地感受到循环的低效并安利一下Mathematica中其它高效的方法。正如wolray的答案中说的，既然选择了使用Mathematica就应该多利用些MMA独有的美妙函数，毕竟如果只是用循环的话C和Fortran之类的语言效率比MMA不知高到哪里去了。。。。。。
当
然我也不是让大家就不用循环了，毕竟很多时候循环的直观性和易读性带来的便利远远比那点效率重要。只是希望大家在循环之前能稍稍想一下，自己的目的是不是
一定要用循环？可不可以用内置函数代替循环？就像上面的几个例子，将循环换成内置函数程序的简洁性和效率都大幅提高，长此以往相信你一定会爱上MMA的~

题外话——关于用编译提速循环
在
MMA中如果一定要使用循环又对效率有一定要求的话，可以选择使用编译，效率能有极大的提高。比如上面的第4个例子使用Complie编译过后的Do循环
用时只有1.86秒，速度提升了将近100倍。如果电脑中有C编译器的话还可以在Compile中加入CompilationTarget ->
"C"选项，速度还能有所提升。编译过后的代码如下：

In[10]:= cf = Compile[{{n, _Integer}, {times, _Integer}},
Mole[{u},
u = ConstantArray[0., {n, n}];
u[[1]] = ConstantArray[1., n];
Do[
Do[u[[i, j]] =
0.25 (u[[i + 1, j]] + u[[i - 1, j]] + u[[i, j + 1]] +
u[[i, j - 1]]),
{i, 2, n - 1}, {j, 2, n - 1}
], {k, 1, times}];
u
]
];
u3 = cf[100, 3000]; // AbsoluteTiming
ArrayPlot[u3, DataReversed -> True, ColorFunction -> "TemperatureMap"]

Out[11]= {1.86055, Null}

前3个例子也都可以通过编译提速很多，这里就不放代码了，有兴趣的同学可以自己动手试一试，如果遇到问题欢迎在评论中与我交流。
需要注意的是编译有很多注意事项，这里推荐一篇写的很好的教程，编译中常见的问题里面都有很好的讲解：怎样编译（Compile）/编译的通用规则/学会这6条，你也会编译
但是一般来讲编译很麻烦，而且再怎么编译效率也很难赶上直接用C，所以个人并不特别建议MMA初学者学习编译。