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

    按：在介紹工具之前先對理論基礎進行必要的回顧是很必要的。沒有理論的基礎，講再多的應用都是空中樓閣。本文主要設涉及線性代數和矩陣論的基本內容。先回顧這部分理論基礎，然後給出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的幫助資訊