科學計算：Python VS. MATLAB(3)----線性代數(shù)基礎

    按：在介紹工具之前先對理論基礎進行必要的回顧是很必要的。沒有理論的基礎，講再多的應用都是空中樓閣。本文主要設涉及線性代數(shù)和矩陣論的基本內容。先回顧這部分理論基礎，然后給出MATLAB，繼而給出Python的處理。個人感覺，因為Python是面向對象的，操縱起來會更接近人的正常思維；而MATLAB大多是以函數(shù)實現(xiàn)的，是向對象施加的一個操作。比如，A是一個矩陣，它有一個屬性attr。用Python更可能是A.attr，而用MATLAB更可能是attr(A)。

一、線形代數(shù)理論基礎

    線形代數(shù)（linear algebra）是數(shù)學的一個分支，研究矩陣理論、向量空間、線性變換和有限維線形方程組等內容。

    比較重要的思想有：1.線性代數(shù)的核心內容是研究有限維線性空間的結構和線性空間的線性變換；2.向量的線性相關性是研究線性空間結構與線性變換理論的基礎；3.矩陣是有限維線性空間的線性變換的表示形式；4.線性方程組的求解問題是n維空間到m維空間線性映射求核和全體原象的問題；5.行列式是研究這些問題的一個工具。

    主要內容有：1.矩陣運算：加減乘除、轉置、逆矩陣、行列式、矩陣的冪、伴隨矩陣；2.矩陣分塊、秩、跡；3.解方程；4.線性相關；5.向量空間；6.特征值和特征向量；7.對稱、相似；8.二次標準型；9.線性空間和基變換；10.正交空間；11.矩陣對角化；13.矩陣分解；14.重要數(shù)字特征。

二、MATLAB的處理

1.建立矩陣

MATLAB中，矩陣是默認的數(shù)據(jù)類型。它把向量看做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內置了很多特殊的矩陣生成函數(shù)，建立特殊矩陣十分方便。

i)第一組用來生成特殊規(guī)則的矩陣。如全零、全一、隨機、等步長等形式。

X=zeros(m,n)   

%生成一個m*n的全0矩陣。同理，ones(m,n)生成一個全1矩陣；eye(m,n)生成一個單位陣。它們的重要作用在于預先分配矩陣空間，所以，在預知矩陣規(guī)模但是不知道矩陣具體數(shù)據(jù)的情況下，先用這幾個函數(shù)生成一個矩陣，對提高運算速度十分有用。

X=rand(m,n)  

%生成一個平均分布的隨機矩陣，數(shù)值區(qū)間[0,1]。同理，randn(m,n)生成一個服從正態(tài)分布的隨機矩陣。注意，這些所謂的隨機實際上都是偽隨機。

v=linspace(a,b,n)   

%產(chǎn)生線性空間矢量。a和b分別是起點和終點，n是本區(qū)間內的點數(shù)，默認100個點。同理，logspace(a,b,n)產(chǎn)生對數(shù)空間矢量。不過它默認點數(shù)是50個。

v=1:0.1:10   

%產(chǎn)生一個線性的矢量。規(guī)格是---起點：步長值：終點

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)  

 %分別產(chǎn)生矩陣X的上三角矩陣和下三角矩陣。

fliplr(X)/flipud(X)/rot90(X)   

%這都是對矩陣的翻轉操作，獲得新的矩陣。分別是左右翻轉（left-right）、上下翻轉（up-down）和逆時針旋轉90°操作。

iii)第三組用來生成一些具有理論價值的，往往是以數(shù)學家命名的矩陣。

magic(n)生成行列相加均為同一個數(shù)字的方陣。pascal(n)生成帕斯卡爾矩陣。hilb(n)生成希爾伯特矩陣。vander(v)生成范德蒙德矩陣。等等。

這些矩陣一般都有相應的學術背景，用到的時候，可以用命令help elmat在最后一個欄目中看看有沒有自己要找的特殊矩陣，如果有，點進去進一步研究即可。

2.矩陣的特征信息

size(X)  %獲得矩陣X的行、列數(shù)。比如，X是一個3*5的矩陣，p=size(X)返回p=[3 5]

length()   %對于矢量，返回的是矢量的長度；對數(shù)組，返回的是數(shù)組最長的那一個維度的長度。

ndims()   %相當于length(size(x))。

numel()   %數(shù)組中元素的個數(shù)。

isempty()和isequal()等is*型函數(shù)    %測試矩陣是否滿足某些條件

[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)   %矩陣或者向量的范數(shù)，具體使用用到再說

c = cond(X)    %矩陣X的條件數(shù)

3.矩陣分解

矩陣分解是矩陣論的重要內容。常用的分解形式在MATLAB中都有函數(shù)予以實現(xiàn)，并且有些分解考慮了多種情況。常見的如：eig()、qr()、schur()、svd()、chol()、lu()等。具體使用的時候

4.矩陣運算

MATLAB默認的是矩陣運算，所以如果想要按元素依次計算，在原來運算符前加一個.號。比如.*表示按元素相乘。

每一個運算符都有一個對應的函數(shù)。如：

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前綴表示的是按照復數(shù)操作進行轉置。

此外，還有一些比較常用的運算：

C=cross(A,B)  

%矢量叉乘。類似的，C=dot(A,B)  是矢量點乘
B = prod(A,dim)   

%數(shù)組元素的乘積，默認按列計算。dim=1是列，dim=2是按行。這個概念很重要！！
類似的，B = sum(A,dim)   求數(shù)組元素的和。dim意義和以上同。
expm()   

%矩陣指數(shù)運算。與此類似的logm(), sqrtm()。其中，funm(A,fun)用來計算矩陣A對通用函數(shù)fun的函數(shù)值。

5.矩陣索引

選擇使用矩陣中的某些元素，就是所謂的矩陣索引了。

A(:,j)   %選取矩陣A的所有行，第j列，同理，A(i,:)是第i行，所有列

A(:,j:k)    %所有行，第j列至第k列（起點和終點均含）

三、Python的處理

Python使用NumPy包完成了對N-維數(shù)組的快速便捷操作。使用這個包，需要導入numpy。SciPy包以NumPy包為基礎，大大的擴展了numpy的能力。為了使用的方便，scipy包在最外層名字空間中包括了所有的numpy內容，因此只要導入了scipy，不必在單獨導入numpy了！但是為了明確哪些是numpy中實現(xiàn)的，哪些是scipy中實現(xiàn)的，本文還是進行了區(qū)分。以下默認已經(jīng)：import numpy as np 以及 impor scipy as sp

下面簡要介紹Python和MATLAB處理數(shù)學問題的幾個不同點。1.MATLAB的基本是矩陣，而numpy的基本類型是多為數(shù)組，把matrix看做是array的子類。2.MATLAB的索引從1開始，而numpy從0開始。

1.建立矩陣

a1=np.array([1,2,3],dtype=int)   #建立一個一維數(shù)組，數(shù)據(jù)類型是int。也可以不指定數(shù)據(jù)類型，使用默認。幾乎所有的數(shù)組建立函數(shù)都可以指定數(shù)據(jù)類型，即dtype的取值。

a2=np.array([[1,2,3],[2,3,4]])   #建立一個二維數(shù)組。此處和MATLAB的二維數(shù)組（矩陣）的建立有很大差別。

同樣，numpy中也有很多內置的特殊矩陣：

b1=np.zeros((2,3))    #生成一個2行3列的全0矩陣。注意，參數(shù)是一個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函數(shù)，即按照某一個已知的數(shù)組的規(guī)模（幾行幾列）建立同樣規(guī)模的特殊數(shù)組。這樣的函數(shù)有zeros_like()、empty_like()、ones_like()，它們的參數(shù)均為如此形式：zeros_like(a,dtype=)，其中，a是一個已知的數(shù)組。

c1=np.arange(2,3,0.1)   #起點，終點，步長值。含起點值，不含終點值。

c2=np.linspace(1,4,10)    #起點，終點，區(qū)間內點數(shù)。起點終點均包括在內。同理，有l(wèi)ogspace()函數(shù)

d1=np.linalg.companion(a)    #伴隨矩陣

d2=np.linalg.triu()/tril()   #作用同MATLAB中的同名函數(shù)

e1=np.random.rand(3,2)    #產(chǎn)生一個3行2列的隨機數(shù)組。同一空間下，有randn()/randint()等多個隨機函數(shù)

fliplr()/flipud()/rot90()    #功能類似MATLAB同名函數(shù)。

xx=np.roll(x,2)   #roll()是循環(huán)移位函數(shù)。此調用表示向右循環(huán)移動2位。

2.數(shù)組的特征信息

先假設已經(jīng)存在一個N維數(shù)組X了，那么可以得到X的一些屬性，這些屬性可以在輸入X和一個.之后，按tab鍵查看提示。這里明顯看到了Python面向對象的特征。

X.flags    #數(shù)組的存儲情況信息。

X.shape    #結果是一個tuple，返回本數(shù)組的行數(shù)、列數(shù)、……

X.ndim   #數(shù)組的維數(shù)，結果是一個數(shù)

X.size    #數(shù)組中元素的數(shù)量

X.itemsize    #數(shù)組中的數(shù)據(jù)項的所占內存空間大小

X.dtype    #數(shù)據(jù)類型

X.T   #如果X是矩陣，發(fā)揮的是X的轉置矩陣

X.trace()    #計算X的跡

np.linalg.det(a)   #返回的是矩陣a的行列式

np.linalg.norm(a,ord=None)    #計算矩陣a的范數(shù)

np.linalg.eig(a)    #矩陣a的特征值和特征向量

np.linalg.cond(a,p=None)    #矩陣a的條件數(shù)

np.linalg.inv(a)    #矩陣a的逆矩陣

3.矩陣分解

常見的矩陣分解函數(shù)，numpy.linalg均已經(jīng)提供。比如cholesky()/qr()/svd()/lu()/schur()等。某些算法為了方便計算或者針對不同的特殊情況，還給出了多種調用形式，以便得到最佳結果。

4.矩陣運算

np.dot(a,b)用來計算數(shù)組的點積；vdot(a,b)專門計算矢量的點積，和dot()的區(qū)別在于對complex數(shù)據(jù)類型的處理不一樣；innner(a,b)用來計算內積；outer(a,b)計算外積。

專門處理矩陣的數(shù)學函數(shù)在numpy的子包linalg中定義。比如np.linalg.logm(A)計算矩陣A的對數(shù)?？梢?，這個處理和MATLAB是類似的，使用一個m后綴表示是矩陣的運算。在這個空間內可以使用的有cosm()/sinm()/signm()/sqrtm()等。其中常規(guī)exp()對應有三種矩陣形式：expm()使用Pade近似算法、expm2()使用特征值分析算法、expm3()使用泰勒級數(shù)算法。在numpy中，也有一個計算矩陣的函數(shù)：funm(A,func)。

5.索引

numpy中的數(shù)組索引形式和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屬性被重新賦值，要求就是元素個數(shù)不變。2*5=10

print x[1,3]    #二維數(shù)組索引單個元素，第2行第4列的那個元素

print x[0]   #第一行所有的元素

y=np.arange(35).reshape(5,7)    #reshape()函數(shù)用于改變數(shù)組的維度

print y[1:5:2,::2]    #選擇二維數(shù)組中的某些符合條件的元素

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

作為第一篇正式的介紹技術操作的文章，終于寫完了，很費勁。恰恰就是在這個費勁的過程中，讓我能更好的認識兩者的區(qū)別和聯(lián)系，同時梳理了展開的思路，摸索出了進一步學習的方法。

我們可以看到，MATLAB中實現(xiàn)了的函數(shù)或者功能，在numpy中都有了對應，并且有些實現(xiàn)的更好。當然，這只是冰山一角。如果你不愿意通讀文檔（很枯燥，誰也不愿意干?。┮矐撚欣碛上嘈?，Python有能勝任工作的實現(xiàn)已經(jīng)存在。后面的內容，將不再這樣列出各種函數(shù)和功能，而是以某一個實際問題為核心，進行專題式的研究。至于全方位的了解，請自己查閱文檔。有個經(jīng)驗之談，就是，應該充分的利用文檔中的【see also】功能，依此追蹤下去，必然會獲得關于某主題的全方位的認識。比如，在查閱ones()的時候，MATLAB的【see also】就給出了complex|eye|true|zeros四個鏈接。這就說明，這幾個函數(shù)其實是有關聯(lián)的，點進去進行簡單的學習，找到共性，那么，可能很多人都遇到過的最大的困惑——那么多函數(shù)怎么記住呀？——就已經(jīng)解決了。因為，我們不需要記住所有的函數(shù)，我們只需要記住有那么回事，只需要記住一個類似的函數(shù)，就可以很輕易的在用的時候順藤摸瓜找出需要的函數(shù)。

下面簡單的給出MATLAB和Python的自查自學方法吧！

1.MATLAB

help 函數(shù)名  

%在控制臺給出某函數(shù)或者主題的幫助信息

doc  函數(shù)名 

%在幫助瀏覽器中給出幫助信息，這個界面更友好。在help browser中既有MATLAB整個產(chǎn)品的瀏覽左窗口，也有一個搜索框。同時還有大量存在的超鏈接。就一個感興趣的主題，點下去，全面學習。不過要記住:別分神哦~~點到最后都忘了自己究竟要做什么！

lookfor 關鍵詞   

%這是一個模糊尋找，含有關鍵詞的詞條入口都會給出來

2.Python

help(np.array)    #查看關于array的幫助信息

help(np.add)  #查看關于add的幫助信息

================update 20121229=================

關于python科學計算，隆重推薦sage math，sage的特點和用法，在本博客較新博文中有介紹。