條件數演算法

⑴ 稀疏矩陣的運算

M AT L A B中對滿矩陣的運算和函數同樣可用在稀疏矩陣中.結果是稀疏矩陣還是滿矩陣,
這取決於運算符或者函數及下列的操作數:當函數用一個矩陣作為輸入參數,輸出參數為一個標量或者一個給定大小的向量時,輸出參數的格式總是返回一個滿陣形式,如命令s i z e.
當函數用一個標量或者一個向量作為輸入參數,輸出參數為一個矩陣時,輸出參數的格式也總是返回一個滿矩陣,如命令e y e.還有一些特殊的命令可以得到稀疏矩陣,如命令s p e y e.
對於單參數的其他函數來說,通常返回的結果和參數的形式是一樣的,如d i a g.
對於雙參數的運算或者函數來說,如果兩個參數的形式一樣,那麼也返回同樣形式的結果.在兩個參數形式不一樣的情況下,除非運算的需要,均以滿矩陣的形式給出結果.
兩個矩陣的組和[A B],如果A或B中至少有一個是滿矩陣,則得到的結果就是滿矩陣.
表達式右邊的冒號是要求一個參數的運算符,遵守這些運算規則.
表達式左邊的冒號不改變矩陣的形式. 假設有:
這是一個5×5的單位滿矩陣和相應的稀疏矩陣.
(a) C = 5*B,結果為:
這是一個稀疏矩陣.
(b) D = A + B,給出的結果為:
這是一個滿矩陣.
(c) x = B h,結果為:
這是一個滿向量.
有許多命令可以對非零元素進行操作.
命令集8 9矩陣的非零元素
n n z ( A )求矩陣A中非零元素的個數.它既可求滿矩陣也可求稀疏矩陣.
s p y ( A )畫出稀疏矩陣A中非零元素的分布.也可用在滿矩陣中,在
這種情況下,只給出非零元素的分布.
s p y ( A , c s t r , s i z e )用指定的顏色c s t r(見表1 3 - 1 )和在s i z e規定的范圍內畫出稀疏
矩陣A中非零元素的分布.
n o n z e r o s ( A )按照列的順序找出矩陣A中非零的元素.
s p o n e s ( A )把矩陣A中的非零元素全換為1.
s p a l l o c ( m , n ,產生一個m×n階只有n z m a x個非零元素的稀疏矩陣.這樣可以
n z m a x )有效地減少存儲空間和提高運算速度.
n z m a x ( A )給出為矩陣A中非零元素分配的內存數.不一定和n n z ( A )得
到的數相同;參見s p a r s e或者s p a l l o c.
i s s p a r s e ( A )如果矩陣A是稀疏矩陣,則返回1;否則返回0.
s p f u n ( f c n , A )用A中所有非零元素對函數f c n求值,如果函數不是對稀疏矩
陣定義的,同樣也可以求值.
s p r a n k( A )求稀疏矩陣A的結構秩.對於所有的矩陣來說,都有
s p r a n k ( A）≥r a n k ( A ). 用下面的命令定義稀疏矩陣:
創建一個大矩陣:
Big=kron(A, A)
這個矩陣B i g是什麼樣子呢?
K r o n e c k e r張量積給出一個大矩陣,它的元素是矩陣A的元素之間可能的乘積.因為參量都是稀疏矩陣,所以得到的矩陣也是一個稀疏矩陣.可以用命令 w h o s和i s s p a r s e來確認一下.
查看矩陣B i g的結構圖,可輸入s p y ( B i g ),結構圖如右圖所示. 從圖中可以看出B i g是一個塊雙對角矩陣. MATLAB中有四個基本稀疏矩陣,它們是單位矩陣,隨機矩陣,對稱隨機矩陣和對角矩陣.
命令集9 0單位稀疏矩陣
s p e y e ( n )生成n×n的單位稀疏矩陣.
s p e y e ( m , n )生成m×n的單位稀疏矩陣.
命令speye(A) 得到的結果和s p a r s e ( e y e ( A ) )是一樣的,但是沒有涉及到滿陣的存儲.
命令集9 1隨機稀疏矩陣
s p r a n d ( A )生成與A有相同結構的隨機稀疏矩陣,且元素服從均勻分布.
s p r a n d ( m , n , d e n s )生成一個m×n的服從均勻分布的隨機稀疏矩陣,有d e n s×m×
n個非零元素,0≤d e n s≤1.參數d e n s是非零元素的分布密度.
s p r a n d ( m , n , d e n s ,生成一個近似的條件數為1 /rc,大小為m×n的隨機稀疏矩
r c )陣.如果rc=rc是一個長度為l≤l ( m i n (m, n) )的向量,那麼
矩陣將rci作為它l個奇異值的第一個,其他的奇異值為0.
s p r a n d n ( A )生成與A有相同結構的隨機稀疏矩陣,且元素服從正態分布.
s p r a n d n ( m , n , d e n s ,生成一個m×n的服從正態分布的隨機稀疏矩陣,和sprand
r c )一樣.
s p r a n d s y m ( S )生成一個隨機對稱稀疏矩陣.它的下三角及主對角線部分與S的結構相同,矩陣元素服從正態分布.
s p r a n d s y m ( n , d e n s )生成一個m×n的隨機對稱稀疏矩陣.矩陣元素服從正態分布,分布密度為d e n s.
s p r a n d s y m ( n , d e n s ,生成一個近似條件數為1 /rc的隨機對稱稀疏矩陣.元素以0r c )對稱分布,但不是正態分布.如果rc=rc是一個向量,則矩陣有特徵值rci.也就是說,如果rc是一個正向量,則矩陣是正定矩陣.
s p r a n d s y m ( n , d e n s ,生成一個正定矩陣.如果k= 1,則矩陣是由一正定對稱矩r c , k )陣經隨機J a c o b i旋轉得到的,其條件數正好等於1 /rc;如果k= 2,則矩陣為外積的換位和,其條件數近似等於1 /rc.
s p r a n d s y m ( S , d e n s ,生成一個與矩陣S結構相同的稀疏矩陣,近似條件數為1 /rc.
r c , 3 )參數d e n s被忽略,但是這個參數在這個位置以便函數能確認最後兩個參數的正確與否. (a) 假設有矩陣A:
輸入R a n d o m = s p r a n d n ( A ),可得到隨機稀疏矩陣:
矩陣中隨機數的位置和矩陣A中非零元素的位置相同.
(b) 對於( a )中的矩陣A,輸入:
B = s p r a n d s y m ( A )
結果為:
這是一個用矩陣A的下三角及主對角線部分創建的對稱矩陣,在非零元素的位置用隨機數作為元素值.
用命令s p d i a g s可以取出對角線元素,並創建帶狀對角矩陣.假設矩陣A的大小為m×n,
在p個對角線上有非零元素.B的大小為m i n (m×n)×p,它的列是矩陣A的對角線.向量d的長度為p,其整型分量給定了A的對角元:
di0 主對角線上的對角線
命令集9 2對角稀疏矩陣
[ B , d ] = s p d i a g s ( A )求出A中所有的對角元,對角元保存在矩陣B中,它們的下標保存在向量d中.
s p d i a g s ( A , d )生成一個矩陣,這個矩陣包含有矩陣A中向量d規定的對角元.
s p d i a g s ( B , d , A )生成矩陣A,用矩陣B中的列替換d定義的對角元.
A = s p d i a g s ( B , d , m , n )用保存在由d定義的B中的對角元創建稀疏矩陣A.
例11 . 4給出了如何使用s p d i a g s命令來解普通微分方程組. 在許多實際應用中要保留稀疏矩陣的結構,但是在計算過程中的中間結果會減弱它的稀疏性,如L U分解.這就會導致增加浮點運算次數和存儲空間.為了避免這種情況發生,在第稀疏矩陣1 2 9
M AT L A B中用命令對矩陣進行重新安排.這些命令都列在下面的命令集9 3中.通過h e l p命令
可以得到每個命令更多的幫助信息,也可見h e l p d e s k.
命令集9 3矩陣變換
c o l m m d ( A )返回一個變換向量,使得矩陣A列的秩為最小.
s y m m m d ( A )返回使對稱矩陣秩為最小的變換.
s y m r c m ( A )矩陣A的C u t h i l l - M c K e e逆變換.矩陣A的非零元素在主對角線附近.
c o l p e r m ( A )返回一個矩陣A的列變換的向量.列按非零元素升序排列.有時這是L U因式分解前有用的變換:lu(A(:, j)).如果A是一個對稱矩陣,對行和列進行排序,這有利於C h o l e s k y分解:chol(A(j, j)).
r a n d p e r m ( n )給出正數1,2,. . .,n的隨機排列,可以用來創建隨機變換矩陣.
d m p e r m ( A )對矩陣A進行D u l m a g e - M e n d e l s o h n分解,輸入help dmperm可得更多信息. 創建一個秩為4的變換矩陣,可輸入:
一旦運行p e r m = r a n d p e r m ( 4 ),就會得到:
給出的變換矩陣為:
如果矩陣A為:
輸入命令:
運行結果為:
有兩個不完全因式分解命令,它們是用來在解大線性方程組前進行預處理的.用h e l p d e s k命令可得更多信息.命令集9 4不完全因式分解c h o l i n c ( A , o p t )進行不完全C h o l e s k y分解,變數o p t取下列值之一:
d r o p t o l指定不完全分解的舍入誤差,0給出完全分解.
m i c h o l如果m i c h o l = 1,則從對角線上抽取出被去掉的元素.
r d i a g用s q r t ( d r o p t o l*n o r m ( X ( : , j ) ) )代替上三角分
解因子中的零元素,j為零元素所在的列.
[ L , U , P ]=返回矩陣X的不完全分解得到的三個矩陣L,U和P,變數o p t取
l u i n c ( X , o p t )下列值之一:
d r o p t o l指定分解的舍入誤差.
m i l u改變分解以便從上三角角分解因子中抽取被去掉的列元素.
u d i a g用d r o p t o l值代替上三角角分解因子中的對角線上的零元素.
t h r e s h中心極限.
解稀疏線性方程組既可用左除運算符解,也可用一些特殊命令來解.
命令集9 5稀疏矩陣和線性方程組
s p p a r m s ( k e y s t r , o p )設置稀疏矩陣演算法的參數,用help spparms可得詳細信息.
s p a u g m e n t ( A , c )根據[ c*l A; A' 0 ]創建稀疏矩陣,這是二次線性方程組的最
小二乘問題.參見7 . 7節.
s y m b f a c t ( A )給出稀疏矩陣的C h o l e s k y和L U因式分解的符號分解因子.
用help symbfact可得詳細信息.
稀疏矩陣的范數計算和普通滿矩陣的范數計算有一個重要的區別.稀疏矩陣的歐幾里德范數不能直接求得.如果稀疏矩陣是一個小矩陣,則用n o r m ( f u l l ( A ) )來計算它的范數;但是對於大矩陣來說,這樣計算是不可能的.然而M AT L A B可以計算出歐幾里德范數的近似值,在計算條件數時也是一樣.
命令集9 6稀疏矩陣的近似歐幾里德范數和條件數
n o r m e s t ( A )計算A的近似歐幾里德范數,相對誤差為1 0-6.
n o r m e s t ( A , t o l )計算A的近似歐幾里德范數,設置相對誤差t o l,而不用預設時的1 0-6.
[ n r m , n i t ] =計算近似n r m范數,還給出計算范數迭代的次數n i t.
n o r m e s t ( A )
c o n d e s t ( A )求矩陣A條件數的1 -范數中的下界估計值.
[ c , v ]=求矩陣A的1 -范數中條件數的下界估計值c和向量v,使得
c o n d e s t ( A , t r )| |Av| | = ( | |A| | | |v| | ) / c.如果給定t r,則給出計算的過程.t r= 1,
給出每步過程;t r=-1,給出商c / r c o n d ( A ). 假設給出:
用n o r m A p p r o x = n o r m e s t ( S p r s )計算出:
用t h e N o r m = n o r m ( f u l l ( S p r s ) )得:
為了找到它們之間的差別,計算d i f f e r e n c e = t h e N o r m - n o r m A p p r o x,得:
在許多應用中,n o r m e s t計算得到的近似值是一個很好的近似歐幾里德范數,它的計算步數要比n o r m要少得多;可參見7 . 6節.
用e t r e e命令來找到稀疏對稱矩陣的消元樹,用向量f來描述消元樹,還可用e t r e e p l o t命令畫出來.元素fi是矩陣的上三角C h o l e s k y分解因子中i行上第1非零元素的列下標.如果有非零元素,則fi= 0.消元樹可以這樣來建立:
節點i是fi的孩子,或者如果fi= 0,則節點i是樹的根節點.
命令集9 7矩陣的消元樹
e t r e e ( A )求A的消元樹向量f,這個命令有可選參數;輸入h e l p
e t r e e獲取幫助.
e t r e e p l o t ( A )畫出向量f定義的消元樹圖形.
t r e e p l o t ( p , c , d )畫出指針向量p的樹圖形,參數c和d分別指定節點的顏色和分支數.e t r e e p l o t可以調用這個命令.
t r e e l a y o u t顯示樹的結構,t r e e p l o t可以調用這個命令. 假設有對稱稀疏矩陣B:
運行命令b t r e e = e t r e e ( B ),結果為:
開始的數字2不難理解,它是矩陣的第1列上第1個非零元素的行數,它決定了在C h o l e s k y分解因子的第1行第2列處有一個非零元素.當縮減第1列的元素時就得到第2列的數字5.B在縮減後,在( 5 , 2 )位置的元素是非零的,這樣消元樹向量中第2個元素的值為5.
s p y ( c h o l ( B ) )給出了C h o l e s k y分解因子的結構圖,如圖9 - 2所示:
圖9-2 Cholesky分解結構圖
圖9-3 矩陣B的消元樹
這個向量消元樹可以這樣來建立:上三角中只有一行有非零元素,節點8,因此這就是樹
唯一的根.節點1是節點2的孩子,節點2和3是節點5的孩子,而節點5是節點6的孩子.節點4和6是節點7的孩子,而節點7又是節點8的孩子,即根的孩子.
命令e t r e e p l o t ( B )給出了樹的結構圖,如圖9 - 3所示.
消元樹的形狀取決於列和行序,它可以用來分析消元過程.
用g p l o t命令可以畫出坐標和矩陣元素間的聯系圖形.必須在n×2的矩陣中給出n個坐標,
矩陣的每一行作為一個點.這樣就創建出點點之間連接的n×n矩陣,如果點4連接到點8,則(4, 8)的值為1.由於是一個大矩陣,而且非零元素較少,所以它應該被建成稀疏矩陣.
這個圖可以說明網路問題,如傳遞問題.它還包含有線性方程組中未知量之間的相關信息.
命令集9 8網路圖形
g p l o t ( A , K )如果矩陣A的a(i, j)不為0,則將點ki連接到點kj.K是一個n×
2的坐標矩陣,A是一個n×n的關聯矩陣.
g p l o t ( A , K , s t r )用字元串s t r給定的顏色和線型畫出的同上圖形.字元串s t r的
取值參見表1 3 - 1.
[ X , A ] = u n m e s h ( E )求網格邊界矩陣E的L a p l a c e矩陣A和網格點的坐標矩陣X.
例七
假設有下面的坐標矩陣K和關聯矩陣A:
矩陣A在稀疏化後,用命令g p l o t ( A , K )畫出圖9 - 4,給出了點(0, 1)和點(4, 1)之間所有可能的路徑.