科學計算：Python VS. MATLAB

儀同3司的博客
http://blog.sina.com.cn/xmatlab
科學計算
http://blog.sina.com.cn/s/articlelist_1596149063_14_1.html
科學計算：Python VS. MATLAB (1)----給我一個理由先
http://blog.sina.com.cn/s/blog_5f234d4701012h9n.html
科學計算：Python VS. MATLAB (2)----準備與前提
http://blog.sina.com.cn/s/blog_5f234d4701012kng.html
科學計算：Python VS. MATLAB(3)----線性代數基礎
http://blog.sina.com.cn/s/blog_5f234d4701012p64.html

按：在介紹工具之前先對理論基礎進行必要的回顧是很必要的。沒有理論的基礎，講再多的應用都是空中樓閣。本文主要設涉及線性代數和矩陣論的基本內容。先回顧這部分理論基礎，然後給出MATLAB，繼而給出Python的處理。個人感覺，因為Python是面向對象的，操縱起來會更接近人的正常思維；而MATLAB大多是以函數實現的，是向對象施加的一個操作。比如，A是一個矩陣，它有一個屬性attr。用Python更可能是A.attr，而用MATLAB更可能是attr(A)。
一、線形代數理論基礎
線形代數（linear algebra）是數學的一個分支，研究矩陣理論、向量空間、線性變換和有限維線形方程組等內容。
比較重要的思想有：1.線性代數的核心內容是研究有限維線性空間的結構和線性空間的線性變換；2.向量的線性相關性是研究線性空間結構與線性變換理論的基礎；3.矩陣是有限維線性空間的線性變換的表示形式；4.線性方程組的求解問題是n維空間到m維空間線性映射求核和全體原象的問題；5.行列式是研究這些問題的一個工具。
主要內容有：1.矩陣運算：加減乘除、轉置、逆矩陣、行列式、矩陣的冪、伴隨矩陣；2.矩陣分塊、秩、跡；3.解方程；4.線性相關；5.向量空間；6.特征值和特征向量；7.對稱、相似；8.二次標準型；9.線性空間和基變換；10.正交空間；11.矩陣對角化；13.矩陣分解；14.重要數字特征。
二、MATLAB的處理
1.建立矩陣
MATLAB中，矩陣是默認的數據類型。它把向量看做1×N或者N×1的矩陣。
%建立了一個行向量，不同元素之間使用空格或者逗號分開都是可以的。
A=[1,2,3] 或者 A=[1 2 3]
%建立一個矩陣，使用分號隔開不同的行。
A=[1,2,3;4,5,6]
%那麽，建立一個列向量就好辦了。每行一個元素，分號分開即可。當然也可以使用行向量的轉置（一個撇號表示轉置）。
A=[1;2;3] 或者 A=[1,2,3]’
MATLAB內置了很多特殊的矩陣生成函數，建立特殊矩陣十分方便。
i)第一組用來生成特殊規則的矩陣。如全零、全一、隨機、等步長等形式。
X=zeros(m,n)
%生成一個m*n的全0矩陣。同理，ones(m,n)生成一個全1矩陣；eye(m,n)生成一個單位陣。它們的重要作用在於預先分配矩陣空間，所以，在預知矩陣規模但是不知道矩陣具體數據的情況下，先用這幾個函數生成一個矩陣，對提高運算速度十分有用。
X=rand(m,n)
%生成一個平均分布的隨機矩陣，數值區間[0,1]。同理，randn(m,n)生成一個服從正態分布的隨機矩陣。註意，這些所謂的隨機實際上都是偽隨機。
v=linspace(a,b,n)
%產生線性空間矢量。a和b分別是起點和終點，n是本區間內的點數，默認100個點。同理，logspace(a,b,n)產生對數空間矢量。不過它默認點數是50個。
v=1:0.1:10
%產生一個線性的矢量。規格是---起點：步長值：終點
ii)第二組用來在原有矩陣基礎上獲得一個具有某些特征的矩陣。
X=diag(v,k)和v=diag(X,k)
%前者用矢量v中的元素生成一個對角矩陣，k是對角移位因子，默認為0，即主對角。k>0，對角線右移。後者返回矩陣X的對角元素，存在矢量v中。k的意義相同。
X1=triu(X,k)和X1=tril(X,k)
%分別產生矩陣X的上三角矩陣和下三角矩陣。
fliplr(X)/flipud(X)/rot90(X)
%這都是對矩陣的翻轉操作，獲得新的矩陣。分別是左右翻轉（left-right）、上下翻轉（up-down）和逆時針旋轉90°操作。
iii)第三組用來生成一些具有理論價值的，往往是以數學家命名的矩陣。
magic(n)生成行列相加均為同一個數字的方陣。pascal(n)生成帕斯卡爾矩陣。hilb(n)生成希爾伯特矩陣。vander(v)生成範德蒙德矩陣。等等。
這些矩陣一般都有相應的學術背景，用到的時候，可以用命令help elmat在最後一個欄目中看看有沒有自己要找的特殊矩陣，如果有，點進去進一步研究即可。
2.矩陣的特征信息
size(X) %獲得矩陣X的行、列數。比如，X是一個3*5的矩陣，p=size(X)返回p=[3 5]
length() %對於矢量，返回的是矢量的長度；對數組，返回的是數組最長的那一個維度的長度。
ndims() %相當於length(size(x))。
numel() %數組中元素的個數。
isempty()和isequal()等is*型函數 %測試矩陣是否滿足某些條件
[V,D] = eig(A) %矩陣A的特征值D和特征向量V。
k = rank(A) %矩陣A的秩
b = trace(A) %矩陣A的跡，即對角線元素之和
d = det(X) %方陣A的行列式
Y = inv(X) %矩陣X的逆矩陣
n = norm(X,option) %矩陣或者向量的範數，具體使用用到再說
c = cond(X) %矩陣X的條件數
3.矩陣分解
矩陣分解是矩陣論的重要內容。常用的分解形式在MATLAB中都有函數予以實現，並且有些分解考慮了多種情況。常見的如：eig()、qr()、schur()、svd()、chol()、lu()等。具體使用的時候
4.矩陣運算
MATLAB默認的是矩陣運算，所以如果想要按元素依次計算，在原來運算符前加一個.號。比如.*表示按元素相乘。
每一個運算符都有一個對應的函數。如：
A+B=plus(A,B)、A-B=minus(A,B)
A*B=mtimes(A,B)、A.*B=times(A,B)
A/B=mrdivide(A,B)、A./B=rdivide(A,B)、A\B=mldivide(A,B)、A.\B=ldivide(A,B)
A^B=mpower(A,B)、A.^B=power(A,B)
A'=ctranspose(A)、A.'=transpose(A)
其中的前綴m自然是表示matrix的意思。沒有m前綴的就是按元素進行的意思。最後那個轉置操作，c前綴表示的是按照復數操作進行轉置。
此外，還有一些比較常用的運算：
C=cross(A,B)
%矢量叉乘。類似的，C=dot(A,B) 是矢量點乘
B = prod(A,dim)
%數組元素的乘積，默認按列計算。dim=1是列，dim=2是按行。這個概念很重要！！
類似的，B = sum(A,dim) 求數組元素的和。dim意義和以上同。
expm()
%矩陣指數運算。與此類似的logm(), sqrtm()。其中，funm(A,fun)用來計算矩陣A對通用函數fun的函數值。
5.矩陣索引
選擇使用矩陣中的某些元素，就是所謂的矩陣索引了。
A(:,j) %選取矩陣A的所有行，第j列，同理，A(i,:)是第i行，所有列
A(:,j:k) %所有行，第j列至第k列（起點和終點均含）
三、Python的處理
Python使用NumPy包完成了對N-維數組的快速便捷操作。使用這個包，需要導入numpy。SciPy包以NumPy包為基礎，大大的擴展了numpy的能力。為了使用的方便，scipy包在最外層名字空間中包括了所有的numpy內容，因此只要導入了scipy，不必在單獨導入numpy了！但是為了明確哪些是numpy中實現的，哪些是scipy中實現的，本文還是進行了區分。以下默認已經：import numpy as np 以及 impor scipy as sp
下面簡要介紹Python和MATLAB處理數學問題的幾個不同點。1.MATLAB的基本是矩陣，而numpy的基本類型是多為數組，把matrix看做是array的子類。2.MATLAB的索引從1開始，而numpy從0開始。
1.建立矩陣
a1=np.array([1,2,3],dtype=int) #建立一個一維數組，數據類型是int。也可以不指定數據類型，使用默認。幾乎所有的數組建立函數都可以指定數據類型，即dtype的取值。
a2=np.array([[1,2,3],[2,3,4]]) #建立一個二維數組。此處和MATLAB的二維數組（矩陣）的建立有很大差別。
同樣，numpy中也有很多內置的特殊矩陣：
b1=np.zeros((2,3)) #生成一個2行3列的全0矩陣。註意，參數是一個tuple：(2,3)，所以有兩個括號。完整的形式為：zeros(shape,dtype=)。相同的結構，有ones()建立全1矩陣。empty()建立一個空矩陣，使用內存中的隨機值來填充這個矩陣。
b2=identity(n) #建立n*n的單位陣，這只能是一個方陣。
b3=eye(N,M=None,k=0) #建立一個對角線是1其余值為0的矩陣，用k指定對角線的位置。M默認None。
此外，numpy中還提供了幾個like函數，即按照某一個已知的數組的規模（幾行幾列）建立同樣規模的特殊數組。這樣的函數有zeros_like()、empty_like()、ones_like()，它們的參數均為如此形式：zeros_like(a,dtype=)，其中，a是一個已知的數組。
c1=np.arange(2,3,0.1) #起點，終點，步長值。含起點值，不含終點值。
c2=np.linspace(1,4,10) #起點，終點，區間內點數。起點終點均包括在內。同理，有logspace()函數
d1=np.linalg.companion(a) #伴隨矩陣
d2=np.linalg.triu()/tril() #作用同MATLAB中的同名函數
e1=np.random.rand(3,2) #產生一個3行2列的隨機數組。同一空間下，有randn()/randint()等多個隨機函數
fliplr()/flipud()/rot90() #功能類似MATLAB同名函數。
xx=np.roll(x,2) #roll()是循環移位函數。此調用表示向右循環移動2位。
2.數組的特征信息
先假設已經存在一個N維數組X了，那麽可以得到X的一些屬性，這些屬性可以在輸入X和一個.之後，按tab鍵查看提示。這裏明顯看到了Python面向對象的特征。
X.flags #數組的存儲情況信息。
X.shape #結果是一個tuple，返回本數組的行數、列數、……
X.ndim #數組的維數，結果是一個數
X.size #數組中元素的數量
X.itemsize #數組中的數據項的所占內存空間大小
X.dtype #數據類型
X.T #如果X是矩陣，發揮的是X的轉置矩陣
X.trace() #計算X的跡
np.linalg.det(a) #返回的是矩陣a的行列式
np.linalg.norm(a,ord=None) #計算矩陣a的範數
np.linalg.eig(a) #矩陣a的特征值和特征向量
np.linalg.cond(a,p=None) #矩陣a的條件數
np.linalg.inv(a) #矩陣a的逆矩陣
3.矩陣分解
常見的矩陣分解函數，numpy.linalg均已經提供。比如cholesky()/qr()/svd()/lu()/schur()等。某些算法為了方便計算或者針對不同的特殊情況，還給出了多種調用形式，以便得到最佳結果。
4.矩陣運算
np.dot(a,b)用來計算數組的點積；vdot(a,b)專門計算矢量的點積，和dot()的區別在於對complex數據類型的處理不一樣；innner(a,b)用來計算內積；outer(a,b)計算外積。
專門處理矩陣的數學函數在numpy的子包linalg中定義。比如np.linalg.logm(A)計算矩陣A的對數。可見，這個處理和MATLAB是類似的，使用一個m後綴表示是矩陣的運算。在這個空間內可以使用的有cosm()/sinm()/signm()/sqrtm()等。其中常規exp()對應有三種矩陣形式：expm()使用Pade近似算法、expm2()使用特征值分析算法、expm3()使用泰勒級數算法。在numpy中，也有一個計算矩陣的函數：funm(A,func)。
5.索引
numpy中的數組索引形式和Python是一致的。如：
x=np.arange(10)
print x[2] #單個元素，從前往後正向索引。註意下標是從0開始的。
print x[-2] #從後往前索引。最後一個元素的下標是-1
print x[2:5] #多個元素，左閉右開，默認步長值是1
print x[:-7] #多個元素，從後向前，制定了結束的位置，使用默認步長值
print x[1:7:2] #指定步長值
x.shape=(2,5) #x的shape屬性被重新賦值，要求就是元素個數不變。2*5=10
print x[1,3] #二維數組索引單個元素，第2行第4列的那個元素
print x[0] #第一行所有的元素
y=np.arange(35).reshape(5,7) #reshape()函數用於改變數組的維度
print y[1:5:2,::2] #選擇二維數組中的某些符合條件的元素
-------------------------------------------------分割線-------------------------------------------------
作為第一篇正式的介紹技術操作的文章，終於寫完了，很費勁。恰恰就是在這個費勁的過程中，讓我能更好的認識兩者的區別和聯系，同時梳理了展開的思路，摸索出了進一步學習的方法。
我們可以看到，MATLAB中實現了的函數或者功能，在numpy中都有了對應，並且有些實現的更好。當然，這只是冰山一角。如果你不願意通讀文檔（很枯燥，誰也不願意幹！）也應該有理由相信，Python有能勝任工作的實現已經存在。後面的內容，將不再這樣列出各種函數和功能，而是以某一個實際問題為核心，進行專題式的研究。至於全方位的了解，請自己查閱文檔。有個經驗之談，就是，應該充分的利用文檔中的【see also】功能，依此追蹤下去，必然會獲得關於某主題的全方位的認識。比如，在查閱ones()的時候，MATLAB的【see also】就給出了complex|eye|true|zeros四個鏈接。這就說明，這幾個函數其實是有關聯的，點進去進行簡單的學習，找到共性，那麽，可能很多人都遇到過的最大的困惑——那麽多函數怎麽記住呀？——就已經解決了。因為，我們不需要記住所有的函數，我們只需要記住有那麽回事，只需要記住一個類似的函數，就可以很輕易的在用的時候順藤摸瓜找出需要的函數。
下面簡單的給出MATLAB和Python的自查自學方法吧！
1.MATLAB
help 函數名
%在控制臺給出某函數或者主題的幫助信息
doc 函數名
%在幫助瀏覽器中給出幫助信息，這個界面更友好。在help browser中既有MATLAB整個產品的瀏覽左窗口，也有一個搜索框。同時還有大量存在的超鏈接。就一個感興趣的主題，點下去，全面學習。不過要記住:別分神哦~~點到最後都忘了自己究竟要做什麽！
lookfor 關鍵詞
%這是一個模糊尋找，含有關鍵詞的詞條入口都會給出來
2.Python
help(np.array) #查看關於array的幫助信息
help(np.add) #查看關於add的幫助信息