python矩阵运算for循环_如何⽤Python科学计算中的矩阵替代
循环
展开全部
因为在Mathematica中使⽤循环确实是低效的。
32313133353236313431303231363533e78988e69d8331333361313961。。。。。
深层次的原因涉及到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
为了便于计时⽤Module封装⼀下，运⾏时间是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；
random在python中的意思
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倍的效率差距
(图⽚太长了这⾥就不直接显⽰了，链接放在下⾯)
===========================我是结尾的分隔线===============================
这个答案其实从⼀开始就跑题了，还写了这么长的⽬的就在于希望让⼤家切实地感受到循环的低效并安利⼀下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}},
Module[{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初
学者学习编译。
已赞过
已踩过<
你对这个回答的评价是？
评论
收起

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。