矩陣和向量運算概覽

內容:
基本線性代數
高級函數
子矩陣和子向量
速度改進

定義:

A, B, C 是矩陣
u, v, w 是向量
i, j, k 是整數
t, t1, t2 是標量
r, r1, r2 範圍(ranges), 例如, range(0, 3)
s, s1, s2 切分(slices), 例如, slice(0, 1, 3)

基本線性代數

標準運算:加法,減法,與一個標量的乘法


C = A + B; C = A - B; C = -A;
w = u + v; w = u - v; w = -u;
C = t * A; C = A * t; C = A / t;
w = t * u; w = u * t; w = u / t;

運算賦值(computed assignements)


C += A; C -= A; 
w += u; w -= u; 
C *= t; C /= t; 
w *= t; w /= t;

內積,外積以及其它乘積


t = inner_prod(u, v);
C = outer_prod(u, v);
w = prod(A, u); w = prod(u, A); w = prec_prod(A, u); w = prec_prod(u, A);
C = prod(A, B); C = prec_prod(A, B);
w = element_prod(u, v); w = element_div(u, v);
C = element_prod(A, B); C = element_div(A, B);

轉化


w = conj(u); w = real(u); w = imag(u);
C = trans(A); C = conj(A); C = herm(A); C = real(A); C = imag(A);

高級函數

范數(norms)


t = norm_inf(v); i = index_norm_inf(v);
t = norm_1(v);   t = norm_2(v); 
t = norm_inf(A); i = index_norm_inf(A);
t = norm_1(A);   t = norm_frobenius(A); 

乘積(products)


axpy_prod(A, u, w, true);  // w = A * u
axpy_prod(A, u, w, false); // w += A * u
axpy_prod(u, A, w, true);  // w = trans(A) * u
axpy_prod(u, A, w, false); // w += trans(A) * u
axpy_prod(A, B, C, true);  // C = A * B
axpy_prod(A, B, C, false); // C += A * B

注意:函數axpy_prod 的最後一個參數 (bool init)是可選的。當前這個參數的值缺省為true,但在以後可能會發生改變。將init 設置為true 等價於在 axpy_prod之前調用w.clear()。到目前為止,這些函數針對於壓縮矩陣進行了特化,相對於函數 prod 產生了巨大的速度提升。


w = block_prod<matrix_type, 64> (A, u); // w = A * u
w = block_prod<matrix_type, 64> (u, A); // w = trans(A) * u
C = block_prod<matrix_type, 64> (A, B); // w = A * B

注意: 塊的大小(blocksize)可以是任意的整數。然而,總體的速度極大地依賴於塊大小(blocksize),CPU和編譯器的組合。函數block_prod 針對於巨大的密集矩陣而定義。

rank-k 更新


opb_prod(A, B, C, true);  // C = A * B
opb_prod(A, B, C, false); // C += A * B

注意:函數axpy_prod 的最後一個參數 (bool init)是可選的。當前這個參數的值缺省為true,但在以後可能會發生改變。如果A的列數比行數小,那麼這個函數可能給出一定的速度提升,因為這個函數通過內積的和(sum of inner product)來計算 。

子矩陣和子向量

通常策略 使用函數project來訪問子矩陣和子向量:


w = project(u, r);         // 由索引範圍 r 所指定的 u 的子向量
w = project(u, s);         // 由索引切分 s 所指定的 u 的子向量
C = project(A, r1, r2);    // 由兩個索引範圍r1和r2所指定的子矩陣
C = project(A, s1, s2);    // 由兩個索引切分s1和s2所指定的子矩陣
w = row(A, i); w = column(A, j);    // 將一個矩陣的行或列作為一個向量

通常策略 使用函數project來賦值給一個子矩陣和子向量:


project(u, r) = w;         // 賦值給由索引範圍r所指定的u的子向量
project(u, s) = w;         // 賦值給由索引切分s所指定的u的子向量
project(A, r1, r2) = C;    // 賦值給由兩個索引範圍r1和r2所指定的u的子矩陣
project(A, s1, s2) = C;    // 賦值給由兩個索引範圍s1和s2所指定的u的子矩陣
row(A, i) = w; column(A, j) = w;    // 將一個矩陣的行或列作為一個向量

注意: 範圍r = range(start, stop) 包含所有的索引istart <= i < stop。一個切分(slice)是更一般的概念。切分(slice)s = slice(start, stride, size)包含索引start, start+stride, ..., start+(size-1)*stride。 stride 可以是0或負數!如果對於一個範圍start >= stop 或對於一個切分(slice) size == 0 ,那麼就不包含任何元素。

向量和矩陣的子範圍和子切分可以使用函數subrangesublice 來直接生成:


w = subrange(u, 0, 2);         // u的子向量的兩個元素
w = subslice(u, 0, 1, 2);      // u的子向量的兩個元素
C = subrange(A, 0,2, 0,3);     // A的子矩陣的2x3個元素
C = subslice(A, 0,1,2, 0,1,3); // A的子矩陣的2x3個元素
subrange(u, 0, 2) = w;         // 賦值給A的子向量的 2 個元素
subslice(u, 0, 1, 2) = w;      // 賦值給A的子向量的 2 個元素
subrange(A, 0,2, 0,3) = C;     // 賦值給A的子矩陣的2x3個元素
subrange(A, 0,1,2, 0,1,3) = C; // 賦值給A的子矩陣的2x3個元素

有很多方式來將一個矩陣的元素作為一個向量來訪問:

matrix_vector_range<matrix_type> (A, r1, r2);
matrix_vector_slice<matrix_type> (A, s1, s2);

注意: 矩陣策略帶有一個矩陣元素的一個序列,並允許你將其作為一個向量來訪問。尤其是matrix_vector_slice 可以使用一個很通用的方式來完成這一點。matrix_vector_range 的作用稍小一些,因為元素必須沿著對角線分佈。

例子: 為了訪問一個矩陣的子列(sub column)的最開始的兩個元素,我們訪問一個步進(stride)為1的行和一個步進(stride)為0的列:
matrix_vector_slice<matrix_type> (A, slice(0,1,2), slice(0,0,2));

速度改進

矩陣/向量賦值

如果你確定左邊的表達式和右邊的表達式沒有使用相同的內存空間,那麼賦值就沒有aliasing。在這種情況下可以指定一個更高效的賦值方式 :

noalias(C) = prod(A, B);

這種方法避免在通常的賦值中所需要生成的一個臨時矩陣。 'noalias' 賦值要求左邊和右邊是大小相似的(size conformant)。

元素訪問

矩陣元素訪問函數A(i1,i2) 或等價的向量元素訪問函數 (v(i) 或 v[i]) ,在應用於一個稀疏矩陣或一個稀疏向量時通常生成'稀疏元素代理(sparse element proxies)'。proxies 訪問元素而不需要擔心令人討厭的(nasty) C++ 問題:引用是 invalidated。

當應用於一個常量對象的時候,這些'稀疏元素代理(sparse element proxies)' 可以更高效地實現。不幸的是在C++裡面沒有辦法區分在一個賦值的左邊還是右邊訪問一個元素。絕大多數情況下,賦值的右邊的元素訪問不會改變元素,因此使用常量代理(const proxies)可能更好一些。 我們可以在訪問矩陣和向量的元素之前使矩陣和向量為常量來做到這一點,例如:

value = const_cast<const VEC>(v)[i];   // VEC 是 V 的類型

出於相同的理由,如果需要訪問多於一個元素,應當使用const_iterator來代替iterator 。可以通過定義配置宏BOOST_UBLAS_NO_ELEMENT_PROXIES來完全關閉更具挑戰性的(daring) '稀疏元素代理(sparse element proxies)'。

控制嵌套乘積的複雜度

進行下面的運算的複雜度 (加法和乘法運算的數量)?

 R = prod(A, prod(B,C)); 

首先複雜度依賴於矩陣大小。因為函數 prod 是可傳遞的(transitive) (不是可交換的),括號順序(bracket order)影響複雜度。

uBLAS 計算表達式而沒有矩陣和向量臨時變量並注意括號結構(bracketing structure)。然而,對於嵌套的乘積避免臨時變量會不可避免地增加複雜度。反之,通過明確地使用使用變量,嵌套乘積的複雜度可以降低。

針對於這個目的,uBLAS 提供了 3 種可靠的方式 :

 temp_type T = prod(B,C); R = prod(A,T);   // 如果T是預先分配的,那麼這就是可選的方法
 prod(A, temp_type(prod(B,C));
 prod(A, prod<temp_type>(B,C));

'temp_type' 很重要。給定相同類型的 A,B,C 。比如 matrix<float>,選擇是很簡單的。然而,如果 value_type 是混合的(mixed) (int 與 float 或 double) 或者矩陣類型是混合的(mixed) (稀疏且是稀疏的),最好的解決方案就不是那麼明顯。完全依賴於你自己!依賴於 A 的數值屬性以及 prod(B,C)的結果。


Copyright (©) 2000-2007 Joerg Walter, Mathias Koch, Gunter Winkler, Michael Stevens
Use, modification and distribution are subject to the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt ).