為什麼Mathematica里矩陣索引這麼慢，有沒有什麼辦法迴避或者改進？

假設我要對一個矩陣進行如下操作：在每兩行（抱歉一開始寫成了列）中各取一個元素進行一次運算，將結果不斷累加直到遍歷所有不同行的所有元素組合（同一行元素之間不運算）。Mathematica代碼可以這樣寫：
1.Do循環

m=100;n=10;
mat=RandomReal[1,{m,n}];
a=0;
Do[a=a+Sinh[Sin[mat[[i1,j1]]]*Cos[mat[[i2,j2]]]], {i1,m-1}, {i2,i1+1,m},{j1,n},{j2,n}];//AbsoluteTiming


2.Flatten@Table

a=Total@Flatten@Table[Sinh[Sin[mat[[i1,j1]]]*Cos[mat[[i2,j2]]]],{i1,m-1},{i2,i1+1,m},{j1,n},{j2,n}];//AbsoluteTiming


在我的電腦上跑，用Do循環大概是1.56s(用Module封裝後有微小提升)，Table大概是1.37s。後者稍快些。
但是在MATLAB里，這種索引計算，卻要快上20多倍：

m=100;n=10;
mat=rand(m,n);
a=0;
tic
for i1=1:m-1
for i2=i1+1:m
for j1=1:n
for j2=1:n
a=a+sinh(sin(mat(i1,j1))*cos(mat(i2,j2)));
end
end
end
end
toc
% Elapsed time is 0.050881 seconds.


或者即便是用類似Flatten@Table的思路生成低效率的三角矩陣，結果依然很快

tic
for i1=1:m-1
for i2=i1+1:m
for j1=1:n
for j2=1:n
b(i1,i2,j1,j2)=sinh(sin(mat(i1,j1))*cos(mat(i2,j2)));
end
end
end
end
sum(b(:));
toc
% Elapsed time is 0.159301 seconds.


之前有帖子提到Mathematica由於底層設計的原因循環總是比較慢，為什麼在 Mathematica 中使用循環是低效的？ - 編程。不過根據我的實驗，在矩陣索引時Do循環和Table不僅寫法很像，速度其實也差不多。當然，如果是不涉及矩陣索引的常規計算比如Sum這種基本函數，用內建函數的效率確實遠比寫循環高很多，但這並不是MMA獨有的特色，在其他高級語言如MATLAB(已實驗)和R(據說)中也是一樣。

所以總結一下，簡單的循環在MMA里可直接用內建函數，對於矩陣索引，寫Do循環只比Table稍慢，但都要比MATLAB慢上許多倍，這實在讓人有點難以接受。請問有沒有什麼辦法可以提高MMA矩陣索引的速度，或者是我有什麼地方使用錯誤，或者MMA里已有另一套哲學來處理此類問題？

樓主要求的是矩陣的外積，用 3 x 2 的符號矩陣說明，就是

In[71]:= m = {{c11, c12}, {c21, c22}, {c31, c32}};
In[72]:= Total@
Sinh@Flatten[
Table[Outer[Times, Sin[m[[i]]], Cos[m[[i + 1 ;;]]]], {i, 2}]]

Out[72]= Sinh[Cos[c21] Sin[c11]] + Sinh[Cos[c22] Sin[c11]] +
Sinh[Cos[c31] Sin[c11]] + Sinh[Cos[c32] Sin[c11]] +
Sinh[Cos[c21] Sin[c12]] + Sinh[Cos[c22] Sin[c12]] +
Sinh[Cos[c31] Sin[c12]] + Sinh[Cos[c32] Sin[c12]] +
Sinh[Cos[c31] Sin[c21]] + Sinh[Cos[c32] Sin[c21]] +
Sinh[Cos[c31] Sin[c22]] + Sinh[Cos[c32] Sin[c22]]


用這個方法算，應該可以避免過多的矩陣元操作。

In[224]:= m = 100; n = 10;
mat = RandomReal[1, {m, n}];
Total@Sinh@
Flatten[Table[
Outer[Times, Sin[mat[[i]]], Cos[mat[[i + 1 ;;]]]], {i,
m - 1}]] // AbsoluteTiming
a = 0;
(Do[AddTo[a, Sinh[Sin[mat[[i1, j1]]]*Cos[mat[[i2, j2]]]]], {i1,
m - 1}, {i2, i1 + 1, m}, {j1, n}, {j2, n}]; a) // AbsoluteTiming

Out[228]= {1.26389, 194754.}


這篇文章里有一些優化（加速）Mathematica / Wolfram 語言程序的常用方法：10 Tips for Writing Fast Mathematica Code

m = 100;
n = 10;
mat = RandomReal[1, {m, n}];

Module[
{},
m3 = Sinh@TensorProduct[Sin@mat, Cos@mat];
m4 = Array[If[#3 &<= #1, 0, 1] , {m, n, m, n}]; Total[m3*m4, Infinity]] // AbsoluteTiming (* {0.0975016, 198665.} *)

這樣寫不僅好看而且比Map還快.

Total@Flatten@
Table[Sinh[Sin[mat[[i1, j1]]]*Cos[mat[[i2, j2]]]], {i1,
m - 1}, {i2, i1 + 1, m}, {j1, n}, {j2, n}] // AbsoluteTiming

(* {1.71291, 198665.} *)


速度的瓶頸在於你要的是個上三角陣 (i2 &> i1), 如果不這麼要求的話:

Module[{},
m3 = Sinh@TensorProduct[Sin[mat], Cos[mat]];
Total[m3, Infinity]]; // AbsoluteTiming

(* {0.0172109, Null} *)


嗯~改變編程思路大概是可能的，但是那恐怕要在徹底理解代碼意圖的基礎上大改，不過演算法的理解是我的弱項呢……總之快20倍就行了是吧：

m = 100; n = 40;
mat = RandomReal[1, {m, n}];

(* True *)


With[{m=100,n=10},
With[{mat=RandomReal[1,{m,n}]},
Table[Sinh[Sin[mat[[i,j]]] Cos[mat[[k,l]]]] If[k&<=i,0,1], {i,m},{j,n},{k,m},{l,n}]]]~Total~-1//AbsoluteTiming