大家用matlab有遇到過哪些槽點？

01-07

前段時間做數值積分，要積一個形如醬紫的函數
沒錯，在0+是在正負無窮間震蕩的，這是可以證明的。
我先嘗試了MATLAB Numerical Integration (Quadrature) 里的所有函數，一律跪，而且跪得一塌糊塗。

我就只說matlab三個核心戰鬥力好了：
quad：自適應Simpson公式，用遞歸，對震蕩函數求積的效率簡直慘不忍睹，最後甩出 Maximum function count exceeded（默認maxfcnt = 10000），還說可能奇異。
quadl：自適應Gauss-Lobatto公式，同上。BTW，lobatto公式是閉型求積公式，為了能計算端點函數值matlab在處理0+的積點時先加了一個eps(superiorfloat(a,b))。
quadgk：自適應Gauss-Kronrod公式。Kronrod求積公式是很好的演算法，對高度震蕩的函數的定積分的數值求積有非常漂亮的結果，Fortran QUADPACK的QAG採用的就是該公式。matlab使用的7-15Kronrod公式，並且用了高度線性化的代碼，而且用的while循環而不是遞歸不僅能夠提速而且能夠突破500的最大遞歸層數限制，說實話確實是非常精彩，但是最後也是達到最大的區間細分數，拋出一個WA。
修改"MaxIntervalCount"之類的參數使函數能夠突破MATLAB自以為是的奇異性診斷也是沒有用的，前兩個遞歸紛紛拋出超過最大遞歸層數，第三個out of memory，這是由於細分的區間全部存儲在一個叫subs的矩陣里，而每個循環在最差情況下會使得subs矩陣的size翻一番，不難想像如果震蕩次數太多而達不到估計誤差的閾值，內存很快就爆了。

好吧，心灰意冷，不妨試試wolfram|alpha吧。
出乎意料，W|A飛快給出了正確答案。上wolfram查mathematica使用的數值積分的演算法，完全嚇得不敢動彈，函數庫完爆matlab。。。
我就只想問問TMW，號稱面向科學計算的matlab為何數值積分能力弱爆？
單一的數據結構使得很多數據結構層次的演算法難以實現，線性化操作對中等及以下規模的運算效率良好，但是一旦數據量超過1e8，內存立馬紅線並且慢得我寧願用循環。
這有一個簡單的題目：各個數位都不含9的所有8位數的倒數和是多少？
求教大神如何寫出高效的matlab代碼。
BTW，我用python和euler-maclaurin公式把時間減到了12秒。

第一個問題求各位都不含9的8位數字倒數和。

鏈接也看了，沒看到MATLAB比其他語言慢多少啊，貌似還是比較快的吧。

第二個積分昨天跟我的一個朋友(winner245，經常逛MATLAB中文論壇的朋友應該認識)討論了下才發現他已經用MATLAB解決了(MATLAB的符號引擎Mupad也能解，但是速度比較慢)，帖子鏈接：提示信息 - Simwe模擬論壇(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析論壇---(手機驗證註冊)

高振蕩函數的數值積分

下面是他的程序：

format long

warning off

a = 1e-20;

b = 1;

n = 20;

f = @(x) 1./x;

q = @(x) log(x)./x;

qpm = @(x) (1 - log(x))./x.^2;

d = (a+b)/2+eps((a+b)/2);

x = linspace(a,b,n).";

u = bsxfun(@power, x-d, 0:n-1);

uprime = bsxfun(@times, bsxfun(@power, x-d, -1:n-2), 0:n-1);

qprime = repmat(qpm(x),1,n);

A = uprime+1i*qprime.*u;

alpha = Af(x);

I = real(u(end,:)*alpha*exp(1i*q(b)) - u(1,:)*alpha*exp(1i*q(a)))

以及他的評價：

「最近讀了些高震蕩函數積分的論文，才發現樣條內插是一種有效的求解高震蕩積分的演算法，最早的針對樣條的演算法是 Filon 演算法，採用的是2階樣條，後來該演算法被推廣到了B樣條

我在12L給出的解法是 Levin 演算法，這種演算法不是樣條法，它本質是一種匹配演算法。下面對這個演算法做一個簡要的定性介紹。
Levin 演算法是先將高震蕩積分的被積函數表示成另一個函數的微分，這樣，高震蕩積分可以馬上通過萊布尼茲公式寫出結果。所以，高震蕩積分就轉變為了解微分方程。但我們並不需要求出微分方程的所有解，而只需要一組特解即可。但特解也不好直接求。於是，需要用到近似。這裡的近似是將特解表示成為 n 個線性無關的基函數的線性組合，n 個組合係數待定。將這組特解代入微分方程，並選取積分區間上n個特定的點，從而得到關於n個組合係數的n元線性方程組。線性方程組很容易求得。這個求解過程是讓n個組合係數和取的n個點去匹配微分方程，所以，匹配法因此而得名。很顯然，基函數和n個點的選取至關重要。比如，基函數必須滿足線性無關性。在Levin 82年和96年兩篇論文里，它用的都是不同次數的多項式。n個匹配點是等間隔選取的。在Levin後，出現了一些改進演算法，比如，基函數選擇為 Chebychev多項式
高震蕩積分的水很深，Levin 的兩篇paper發表於82和96年，而直到今天，mathematica最新版的高震蕩積分演算法依然是這兩篇論文的演算法。雖然，最近10年出現了一些新論文，但基本上都是小的改進，沒有大的飛躍。

」

第三個indexing的問題，題主的MATLAB方法和numpy的效率相當，我這裡換了另一種方法來和題主的做個對比：(不過numpy的indexing看著的確友好)

I = 100; J = 200; K = 800; L = 300;

v = rand(I, J, K);

idx = unidrnd(K-L, I, J);

idx_k = bsxfun(@plus, idx, reshape(0:L-1, [1,1,L]));

tic;[idx_i, idx_j, ~] = ndgrid(1:I, 1:J, 1:L);

ind = sub2ind([I, J, K], idx_i, idx_j, idx_k);

r = reshape(v(ind(:)), [I, J, L]);toc

%---------------------分割線---------------------------

tic;ind = bsxfun(@plus,(idx-1)*I*J+bsxfun(@plus,(1:I)",0:I:I*(J-1)),I*J*reshape(0:L-1,1,1,L));

rr = v(ind);toc

isequal(r,rr)

Elapsed time is 0.339555 seconds.

Elapsed time is 0.184700 seconds.

ans =

效率提升了近1倍。

說了那麼多，其實我對MATLAB的吐槽還是有不少的。

首先功能太全了，什麼科學計算、統計分析、圖像處理、遙感分析、計算機視覺、金融分析、汽車控制建模、通信、資料庫、信號處理、生物醫學、神經網路、機器學習，現在最新的MATLAB也支持Hadoop和mapreduce了，大規模的並行計算也早就支持了。而且作為一個科學計算平台，和其他語言介面豐富，和Java完美無縫混合，C#和C/C++互相調用也非常方便,最新版的對python的支持也強大了許多，一般來說搞演算法的基本在MATLAB框架下就夠了，實在不方便如果別的語言有現成庫，拿過來用就行，總之，用的越深，中毒越深，越懶得用其他語言，這點非常不好。

其次，MATLAB在國內的教育太垃圾了。大家都知道大學C/C++的教育和實際工作中的應用要求相差甚多，MATLAB就更是了，很多高校雖然都開了MATLAB課，但能把MATLAB教好的幾乎沒有，退一步說，不求學到實用的本領，學生學完課程能對MATLAB有個客觀正確的認識就不錯了。很多人居然認為MATLAB不能用於生產，只能小打小鬧，做些演算法研究的工作。

最後，MATLAB不適合做頻繁細粒度的計算，這樣的工作還是交給C/C++比較好。中大規模以上的運算內存不是問題，MATLAB的內存機制自動化程度很高，對使用者要求低多了。使用者只需平衡函數overhead，數組對象等的訪問開銷和實際計算的開銷比例就行。數據量很大時可以分塊向量化操作，效率一樣杠杠的。我們項目處理的數據量動不動就在幾個T級別，都是分塊處理的，瓶頸根本不在MATLAB這裡。

我遇到過用for是運行正常, 改成parfor,就運行錯誤的, 原因不明.

槽點啊，我覺得最大的槽點是MATLAB自帶編輯器：為什麼不用Emacs或Vim。

另一個槽點是代碼量，很羅嗦。比如，MATLAB和Mathematica，按我經驗，實現相同功能大概是10：1的代碼量。

再有一個是代碼效率，樓上討論很詳盡，這邊不敘述了。另外一提，不同OS上相同軟體相同配置的運行效率貌似也不太一樣。我沒有具體測試過，但記得有一次同一段Mathematica程序在Windows和Linux下運行時間大概是4：1的樣子。

好吧，寫得太水了，樓上大神太多。。

先說個最近遇到的例子吧，是個簡單的定積分

解析解是算不了的。本來也就沒指望它，數值的呢，算了好幾秒，而且得到的結果精度不太好。

fun = @(x,y)(x.^2 + (y - 1/2).^2 &<= 1/4) .* ... ((x - 1).^2 + (y - 1).^2 &<= 1) .* ... ((x - 1).^2 + y.^2 &<= 1); integral2(fun, 0, 1, 0, 1)

試了下Mathematica，算解析解的時間基本上和Matlab算數值解的時間一樣，算數值解瞬間就出來了，而且精度更高

NIntegrate[Boole[x^2+(y-1/2)^2&<1/4(x-1)^2+(y-1)^2&<1(x-1)^2+y^2&<1], {x,0,1},{y,0,1}]

早一段算的一個例子，求一個三元非線性方程組的數值解，Mathematica的速度是Matlab的幾十倍，FindRoot vs fsolve

function r=eqs(a) x=a(1); y=a(2); z=a(3); r=[(-1)*((-1)+x+y)^(-3)+(-1)*y^2*((-1)+x+y)^(-3)+(-13/10)*((-1)+x+y)^(-2)+2*y*((-1)+x+y)^(-2)+(1/10)*((-1)+x+y)^(-3)*(7+13*x+13*y)*z+(-1)*((-1)+x+y)^(-3)*z^2+(1+3*x+y)*(1+z) ^(-2)+(-23/10)*(1+z)^(-1); (-1)*((-1)+x)*((-1)+x+y)^(-3)*((-2)+2*x+y)+(1/10)*((-1)+x+y)^(-3)*((-3)+13*x+13*y+(-10)*z)*((-1)+z)+(1/5)*(1+z)^(-2)*(3+5*x+5*y+8*z)+(1+z)^(-1)*((-9/10)+y*(1+z)^(-1)); (1/10)*(3+(-13)*x+(-13)*y)*((-1)+x+y)^(-2)+((-1)+x+y)^(-2)*z+(-1)*(5+3*x^2+2*((-1)+y)*y+2*x*(1+y))*(1+z)^(-3)+(1/10)*(76+23*x+(-7)*y)*(1+z)^(-2)]; end

sols=fsolve("eqs", [1 -1 0])

FindRoot[{-(y^2/(x + y - 1)^3) + (2 y)/(x + y - 1)^2 - z^2/(x + y - 1)^3 + ((13 x + 13 y + 7) z)/(10 (x + y - 1)^3) - 23/(10 (z + 1)) - 13/(10 (x + y - 1)^2) + (3 x + y + 1)/(z + 1)^2 - 1/(x + y - 1)^3, -(((x - 1) (2 x + y - 2))/(x + y - 1)^3) + ((13 x + 13 y - 10 z - 3) (z - 1))/(10 (x + y - 1)^3) + (5 x + 5 y + 8 z + 3)/(5 (z + 1)^2) + (y/(z + 1) - 9/10)/(z + 1), (-13 x - 13 y + 3)/(10 (x + y - 1)^2) + z/(x + y - 1)^2 + (23 x - 7 y + 76)/(10 (z + 1)^2) - ( 3 x^2 + 2 (y + 1) x + 2 y (y - 1) + 5)/(z + 1)^3 }, {{x, 0}, {y, 0}, {z, 0}}]

並行計算的時候代碼無法調試。。。

難道我的圖形擬合度不好是因為quadgk的問題？？？

有沒有for循環直接關係到運算速度

你要積分的是什麼函數呢？

第二個問題，各個數位都不含9的所有8位數的倒數和是多少？

這種方法只能算到7位的，8位的我的內存不夠。。

a=0:8; [x1,x2,x3,x4,x5,x6,x7]=ndgrid(1:8,a,a,a,a,a,a); A=1/(1000000*x1 + 100000*x2 + 10000*x3 + 1000*x4 + 100*x5 + 10*x6 + x7); sum(A(:))

另外的方法就是寫8個for循環，其實速度不算慢

最近也是剛剛接觸MATLAB，用quadv算貝塞爾函數的積分，提示矩陣奇異，而且得出來的結果簡直慘不忍睹，看來辛普森積分還是要慎用