幾個可以處理矩陣的計算機: octave/rlab/R

---

簡介

要用自由軟體處理矩陣運算, 有好幾個選擇。

Octave [1] 可能是首選, 因為它的語法幾乎與 matlab 相容, 函式庫豐富, 使用社群活絡。 Octave 的文件以 gnu info 格式寫成, 建議一併安裝 pinfo 套件, 方便下 pinfo octave 命令來看手冊。 Windows 用戶請注意: 根據 這頁的說明 您可選擇 (a) 獨立安裝 這個版本 或 (b) 安裝 { cygwin 或 XLiveCD } 及 這個版本; 但這兩個選項彼此衝突。

RLaB [2] 語法比 matlab 漂亮, 是筆者的最愛; 可惜軟體已不再更新。 (也有好處: 不需要追新) 又, 有人延續 R, 發展 rlabplus

R [3] 主要用於統計運算; 矩陣運算只是其中一小部分的功能。 搜尋此套件時, 用 R-base 當關鍵字。

Giac/Xcas [4] 主要用於符號運算 (symbolic computation); 矩陣運算只是其中一小部分的功能。 它的語法比較不同, 另文 介紹。

這幾套軟體由於程式原始碼公開, 都有 *BSD 版, Linux 版, Windows 版...等等, 使用者可以確保學習投資及辛苦設計的程式不會隨著流行潮的改變而付諸流水。 這些軟體都有現成編譯好的執行檔, 可以上 rpmfind 找到。 四者都提供命令列介面, 可使用方向鍵及 「GNU readline 使用者界面」 所定義的一些快速鍵, 其實比圖形介面方便。 R 與 Giac/Xcas 則還另外提供圖形介面。

純量運算

這些軟體的語法雖然有些不同, 但處理矩陣的觀念是一致的, 所以一併介紹。 以下用不同的顏色表示不同軟體的語法; 沒有特別說明的範例就是四者都適用。 現在不妨就請馬上打開一個終端機視窗, 邊看本文邊照打。

每次執行 rlab 時, 它都會去家目錄下找一個文字檔 .rlab 來執行, 如果沒有這個檔案, 它會抱怨找不到 init script。 所以第一次進入 rlab 之前, 請先建立 .rlab 裡面可以是空的; 或者簡單放一句:

        pi = atan(1)*4;

1. 在 shell 下以 octave 或 rlab2 或 R 或 cas 命令進入後, 你就可以開始把它當一部計算機來用: 試著打 3 * (7 - 2), 然後按 Enter, 它就顯示運算結果。 空格多寡無所謂; 但最好把完整的運算式打在同一列上 (octave 與 R 比較不挑剔, 會給你換列繼續打完的機會)。
2. 運算結果可以存在變數裡 (大小寫有別), 例如 x = 3 * (7 - 2) 等一下再拿出來用: 1/x
3. 可以把好幾句話放在同一列, 用分號分開就行:
   x = 3*(7-2); y = 6.2e-3; z = 4.3/8; [octave: 分號有個副作用: 那一句話的運算結果不會顯示出來。]
4. 除了四則運算以外, 常見的數學函數也都可以用, 例如 sin(pi/3) 與 sqrt(3)/2 都會印出 0.866 。 [rlab: pi 不是內建常數, 要自己算, 例如下 pi=atan2(1,1)*4] 三角函數都是以徑度量為單位。 又如 log(9) 會印出 2.1972, 而 exp(2) 則是 7.3891 因為指數對數都是以 2.71828 為底。
5. 也可以處理複數, 例如 sqrt(3+4i) 其中 i 就是 sqrt(-1) 複數的單位。 但要注意: 虛部係數不可省略 (也就是說單獨一個 i 還是應寫成 1i) 且係數與 i 之間不可有空格。
6. 數學運算當中偶爾會出現 1/0 或 0/0 等等狀況, 像這些處理大量數字的程式, 不適合讓這種狀況打斷程式邏輯。 因此它們提供兩個特殊的值: inf 或 infinity 表示無窮大; nan 或 undef 表示 "not-a-number", 確保運算無論如何都有結果。
7. help 命令可以查手冊, 例如 help cross 在 rlab 中會顯示 cross 函數的原始碼; 在 octave 當中則會顯示短短的使用說明。 如果要看比較完整的手冊, 建議在 shell 底下, 下 pinfo octave 或下 rpm -qil rlab | less 並找到 rlab-ref.html
8. 用 help 命令看手冊時, 你面對的是 less 的操作界面, 搜尋字串時甚至可以使用比較複雜的 regular expression [5]。 在 octave (及某些版本的 rlab 中), 如果運算結果太冗長, 一個螢幕放不下 (例如下面提到的矩陣運算), 也是用 less 來讓你分頁瀏覽。
9. 要離開這些數學工具, 回到 shell 底下, 可以按 ^d 鍵 或下 quit 命令。

使用向量與矩陣

建立矩陣的語法: 方括弧內用分號分隔各列; 用逗號分隔一列內各元素, 例如 A = [1,3,-13,4; 5,-6,0,1; -1,-7,5,2] 所代表的就是下面這個矩陣:

                1          3        -13          4  
                5         -6          0          1  
               -1         -7          5          2  

為方便下面舉例, 我們也建一個列向量 b = [-3,5,-19,7] 一個行向量 c = b' 並且呼叫亂數函數 rand 產生一個大矩陣 (裡面的值呈均勻分佈, uniform distribution) x = rand(8,6) 這裡我們看到單引號代表轉置 (transpose)。

取出	[octave 語法]	[rlab 語法]
A3,2	A(3,2)	A[3;2]
b3	b(3)	b[3]
c4	c(4)	c[4]
x 的第 5 列	x(5,:)	x[5;]
x 的第 2 行	x(:,2)	x[;2]
x 的左下角 1/4 小塊	x(5:8,1:3)	x[5:8;1:3]

[rlab: 對矩陣也可以使用「一個註標」的語法, 不過取到的是以 column-major 方式來數的元素, 例如這裡若下 A[7] 會取到 -13。] 這個表示法不僅可以當做 r-value, 也可以當做 l-value。 也就是說, 你不僅可以從矩陣當中的指定位置讀值出來, 還可以把值填入矩陣的指定位置。

其實冒號本身就可以用來製造列向量, 例如 5:8 效果相當於 [5,6,7,8] 甚且可以指定間隔, 例如要製造出 [50,55,60,65,70,75,80] 可以這麼寫: [octave: 50:5:80] [rlab: 50:80:5]

---

以下正在修改當中...

---

1. 整個矩陣裡面的每個元素都乘以十, 然後無條件捨去, 變成 0 到 9 之間的亂數: x = floor(x*10)
2. Q: 請解釋下句的意義。

           [octave 語法]           [rlab 語法]
           x([12,1:4],[6,5])       x[12,1:4; 6,5]
   

3. Q: 如何取得 x 的所有偶數列, 偶數行元素?
4. 你還可以用小矩陣堆成大矩陣 (當然行數或列數要一致), 例如 y = [x;x] 產生一個 16x6 的矩陣; 而 y = [x,x] 則產生一個 8x12 的矩陣。
5. 程式中如何得知一個矩陣 x 的大小? 可用 size(x) 一次取得列數與行數 (傳回值是一個 1x2 的矩陣)。 如果是 rlab, 還可用 x.nr 與 x.nc 分別直接取得 x 的列數與行數
6. 要把整個矩陣 (例如上面的 a) 的所有 行 按順序直著串起來, 變成一個很高的行向量, 在 octave 中用 a(:) 而在 rlab 中用 a[:] 得到的就是 [1;5;-1;...4;1;2] 這個行向量。
7. 那如果想把所有 列 按順序橫著串起來成為一個很寬的列向量呢? 你可以先轉置, 再把所有行串起來, 最後轉置回來, 像這樣: 在 octave 中用 t=a'; t=t(:)' 在 rlab 中用 t=a'[:]' 事實上這個方法可以應用在很多場合, 因為有很多運算只對行做, 或只對列做。
8. Q: sum(x) 會對 x 的每一行分別加總, 產生一個列向量。 請問如何對 x 的每一列分別加總, 產生一個行向量?
9. Q: 請查手冊, 用 max, min 兩函數, 分別求出 x 的奇數列, 各列最大值/各列最小值。

程式實例: Gauss-Jordan 消去法求 reduced row-echelon form

接下來我們用一個完整的範例程式來介紹 octave 與 rlab 的語法。 我們舉的例子是 Gauss-Jordan Elimination。 如果將它拿來用在線性聯立方程組的 係數-常數項 augmented matrix 上, 所得到的就是方程組的解; 而如果用在一個普通的矩陣上, 所得到的新矩陣就叫做原矩陣的 reduced row-echelon form。 以下說明假設讀者已熟知 Gauss-Jordan Elimination 的步驟; 請讀者視需要複習線性代數課本。 (其實只要你會用 「加減消去法」 解線性聯立方程組也就夠了, 我們這裡的目的是用實際的程式範例教你 octave 與 rlab 的語法; 線性代數的理論細節並不那麼重要.) 閱讀本節時, 請務必對照程式碼。

如何載入/執行 octave 版的 lademo.m? Octave 並不搜尋目前目錄下的 .m 檔。 你可以在進入 octave 之後, 寫出完整的路徑用 source 載入 lademo.m, 像這樣:

        source("/home/ckhung/public_html/b/ma/lademo.m");

或是 (一樣是進入 octave 之後) 先設定搜尋目錄, 再打檔名 (那就不要路徑不要加附檔名了), 像這樣:

        LOADPATH="/home/ckhung/public_html/b/ma:";
        lademo

或者更簡單地寫 LOADPATH=".:" 把目前目錄加入搜尋路徑。 注意後面的冒號不可以省略, 否則它會找不到內建函數如 rand 與 sort 等等。

如何載入/執行 rlab 版的 lademo.r? 你同樣可以用指定搜尋路徑的方式來令 rlab 載入你的程式 (下 rlab2 命令之前, 在 shell 下設定環境變數 RLAB2_PATH); 不過因為 rlab 會在目前目錄下搜尋, 所比較簡單的方法是: 在下 rlab2 之前切換到程式檔的目錄, 或是下 rlab2 後用 cd("/home/ckhung/public_html/b/ma/") 切換目錄; 進入 rlab2 之後則可下 rfile lademo 命令 (不要附檔名) 或呼叫 load("lademo.r") 函數來載入。

不論是 octave 或 rlab, 載入我們自己寫的程式之後, 就可以把裡面定義的函數拿來用了:

    A = [1,3,-5,2,94; 4,-3,1,5,58; 6,-2,2,4,72; 0,2,3,-7,-69]
    T = pivot(A,3,2)

octave 版的 lademo.m 與 rlab 版的 lademo.r 程式架構相同, 以下我們只解釋這個共同的架構; 語法請自行參考程式及手冊, 相信這樣會比列舉 if, for, ... 的文法規則更能搔到癢處。 比較特殊的語法, 下兩節再詳述。

Gauss-Jordan Elimination 的基本動作是 pivoting -- 選取一個元素作為 pivot element, 將整列除以它 (所以它就變成了 1), 然後將這一列 (稱為 pivot row) 乘以適當的常數加到上下各列去, 目的在把與 pivot element 同行的所有元素 (合稱為 pivot column) 都清成 0。 程式中的 pivot 函數接受四個參數: 要對那個矩陣做 pivot, pivot 元素在第幾列, pivot 元素在第幾行, 及希望清除那幾列上的對應元素 (因為有時候我們只希望清除 pivot column 上的部分元素, 而不是要對每一列都做這個動作)。 例如 T = rand(5,6) 用亂數產生 T 之後, 可以下 pivot(T,3,4,[2,5]) 以第三列第四行為 pivot 將第二列及第五列的對應元素清除。

程式中的主要函數為 r_r_echelon, 主要有兩個迴圈。 第一步由左而右逐次尋找新的 pivot element 並呼叫 pivot 以清除矩陣的左下角 (見第一個迴圈內變數 c 的變化範圍)。 當然 pivot 不得為 0, 所以必須在第 c 行內上下搜尋以非 0 數字帶頭的列 (迴圈中的變數 p), 並把找到的列整列調上來。 第二步由下而上尋找 pivot element 以清除矩陣的右上角 (見第二個迴圈內變數 r 的變化範圍)。 同樣地, 此時必須在第 r 列內由左而右找到帶頭 (第一個非 0) 的那一行 (迴圈中的變數 p) 才能作為 pivot element。 這裡的動作看起來有點不對稱, 主要是為了能夠處理任意 rank 的矩陣。 另外, 如果還希望顧及 numerical stability, 程式其實很容易可改成 "選取絕對值最大元素作為 pivot"。

我們可以用亂數來產生測試用的矩陣; 但是除非是遇到鬼, 或是中了樂透頭獎 (意思是一樣的 -- 機率近乎 0), 產生的矩陣幾乎必定是 full rank, 它的 reduced row echelon form 幾乎必定呈「左上方塊單位矩陣」的形式, 也因而很難測試出 r_r_echelon 當中比較複雜的邏輯。 為此我們寫了一個函數 rand_mat, 用來產生大小 m 乘 n, 並且 rank 可由使用者指定的亂數矩陣。 這個函數的邏輯稍複雜, 我們略過。 它會用到 rand_permute 產生 1 到 n 的整數的亂數排列。 Q: 請看註解, 學會用 "一句話" (octave 需要兩句話) 呼叫 rand_permute 來幫你做電腦選號。 不過相信數學的人不玩樂透; 而相信自己運氣的人一年不玩超過三次的樂透, 所以請不要問我這個程式準不準 -- 當然是五百廿四萬分之一嘍 :-)

另外還有一個 simplex 函數用來解線性規劃的問題; 以及 lademo.m 當中的 mini 與 lademo.r 當中的 smooth 等輔助函數我們就略過不談了, 請有興趣的讀者自行研究。 我們直接跳到最下面看一小段註解起來的測試程式。 測試時請自己把井號拿掉。 為了方便除錯, 第一句話以一個固定的值當做亂數的 "種子", 這樣每次產生的亂數序列才會相同; 至於 "種子" 的值是多少並不重要。 除錯完畢之後, 產生種子的這句話其實可以刪掉。 然後用亂數產生一個 5x6 的矩陣 (它的秩只有 4), 試一下以第三列第四行為 pivot, 消去第二列及第五列的對應元素, 最後將原矩陣化為 reduced row echelon form。

當然如果不是為了學 octave/rlab 或是學線性代數, 只是想要求方陣 A 的反矩陣, 或是要解聯立方程組 A * x = b 當中的未知行向量 x, 其實只要下 A\b 就可以了, 根本不必如此大費周章 -- 要站在巨人的肩膀上, 不要重新發明輪子, 不是嗎? 上面這個矩陣除法的語法看起來有點奇怪; 不過想想也有道理, 因為矩陣乘法不可交換, 所以除法有分左除右除。 如果 A * x = b, 那麼 x 不就應該是:

                b
              -----
              A

如果你還是不習慣, 在 rlab 當中還可以這樣下: solve(A,b)

Octave 的怪癖與特異功能

1. 你自建的函式庫通常是以「一檔案一函數」的方式來儲存, 且檔名要與函數名稱相同。 我偏好將很多相關函數放入同一個檔案內, 這時候檔內第一句話就不可以 function 開頭, 例如 lademo.m 內以一句 (沒什麼意義的) 1; 開頭
2. Octave 提供很多全域變數讓你更改內定的行為, 例如發生錯誤時要不要「嗶」一聲? 可以設定 beep_on_error 的值來決定。 詳見手冊的 "Variables" 章的 "Summary of Built-in Variables" 篇。
3. Octave 盡力做到與知名的 matlab 相容。 如果你想要從 matlab 跳脫出來, 進入自由軟體的世界, 但是又不希望把過去的 matlab 程式丟掉, 可以參考 FAQ 裡面的 "Porting programs from MATLAB to Octave" [6]。
4. 除了內建函數及隨著 octave 一起安裝的許多工具箱 (見 rpm -qil octave | grep '\.m$' 或手冊) 之外, 還可以在 matlinks [7] 上找到許多移植自 matlab 的數學工具箱。
5. 可以將矩陣存在一個純文字檔裡面 (例如叫做 abc.txt) 內容像這樣:

            2  1  5 -3
           -1  0  4  7
            3 -5  0 -2
   

   然後在 octave 內下 load "abc.txt" 就會將資料讀入 abc 這個矩陣。 附檔名叫什麼其實都無所謂。 這個功能最適合與 awk 或 perl 等文字處理工具配合使用 -- 例如想用 octave 分析其他軟體產生的資料, 可以先用 awk 或 perl 處理成如上格式。

RLaB 的怪癖與特異功能

1. 建議生手從 help INTRO 開始看起。 另外, 下 rpm -ql rlab | grep html 找出手冊來研究, 還可以學到很多進階功能 (例如透過 pipe 與其他程式溝通)
2. rlab 的函數可以接受不定個數的參數, 像是範例中的 pivot 函數的最後一個參數就是。 呼叫者 (caller) 究竟有無提供這些可省略的參數, 在被呼叫者 (callee) 當中, 可以用 exist 函數來檢查。 Octave 雖然也可以傳遞不定個數的參數, 但使用起來麻煩多了。
3. 上述 octave 從純文字檔讀矩陣的功能, rlab 也有。 資料檔的格式一模一樣; 下的則是 x = readm("abc.txt") 命令, 變數的名稱不必與檔案名稱相同。 如果要再讀一次, 必須先用 close 關閉檔案, 否則會讀到空矩陣。
4. rlab 內部計算具有相當高的精準度; 但是印出來的只有 3 位有效數字。 欲將列印時的精準度增加至 7 位, 可下 format(7) 命令。
5. 在 rlab 裡面, 函數的參數不僅可以是數字/字串的矩陣, 還可以是另外一個函數! 例如下面 integrate 函數以 Simpson's method 求積分:

           # 如果看不懂, 可以把每列後面的分號拿掉, 計算過程就會印出來。
           integrate = function(f, low, high, n) {
               if (! exist(n)) { n = 7; }
               h = (high - low) / (2 * n);
               x = (0 : 2*n) * h + low;
               k = (mod( (0 : 2*n), 2*ones(1,2*n+1) ) + 1) * 2;
               k[1,2*n+1] = [1,1];
               return sum(k .* f(x)) * h / 3;
           }
   

   可以這麼用: integrate(exp, 1, 3)

繪圖

Octave 與 rlab 都可以透過 gnuplot [8] 來將向量或矩陣裡面的資料畫出來。 例如要畫 y 相對於 x 的變化趨勢, 可以呼叫 plot 函數:

        [octave 語法]           [rlab 語法]
        x = -pi*5:0.3:pi*5;     x = -pi*5:pi*5:0.3;
        y = sin(x)./x;          y = sin(x)./x;
        plot(x, y)              plot(x, y)

當然 x 與 y 可以不必是解析函數, 也可以是資料。 例如若有一個資料檔 t.txt 內容如下:

   -2.00000   10.93063
   -1.00000    6.12849
    0.00000    0.98502
    1.00000   -0.36626
    2.00000    1.49372
    3.00000    2.55317
    4.00000   11.51355
    5.00000   18.20215
    6.00000   28.11316

我們先讀入原始資料, 並畫出來:

        [octave 語法]                   [rlab 語法]
        load "t.txt"                    t = readm("t.txt")
        x = t(:,1); y = t(:,2);         x = t[;1]; y = t[;2];
        plot(x,y)                       plot(x,y)

接著我們分別試著找一條直線 y = p₁ + p₂ x 及一條二次曲線 y = q₁ + q₂ x + q₂ x₃ 來近似原始資料。 你可以看出用 octave 與 rlab 處理矩陣相關問題有多麼方便: 三言兩語就解決了兩題 normal equations A^T * (A p - y) = 0 與 A^T * (A q - y) = 0 :

        [octave 語法]                   [rlab 語法]
        A = [ones(9,1), x]              A = [ones(9,1), x]
        p = (A'*A) \ A'*y               p = (A'*A) \ A'*y
        A = [ones(9,1), x, x.*x]        A = [ones(9,1), x, x.*x]
        q = (A'*A) \ A'*y               q = (A'*A) \ A'*y

如果你忘記最小方差問題 (least squares) 的解法也沒有關係, 這個例子的重點是如何在一張圖裡面同時呈現三個函數:

        # octave 接受兩個參數, 第一個是三函數的共同 x 座標
        plot(x, [y, p(1)+p(2)*x, q(1)+q(2)*x+q(3)*x.*x])
        # rlab 則接受一個大矩陣作為參數, 第一行是三函數的共同 x 座標
        plot([x, y, p[1]+p[2]*x, q[1]+q[2]*x+q[3]*x.*x])

順便一提, 這些資料是在 octave 裡面用常態分佈亂數產生的:

        x = -2:6
        y = 1.2 * x .* x - x * 2.5 + 0.7
        t = [x; y+normal_rnd(0,0.7,1,9)]'
        save -ascii "data.txt" t

繪圖的功能不限於平面圖形。 3-d 網格, 不論是函數或資料點, 也都可以畫。 例如想畫一個馬鞍曲面 z = x*x - y*y 可以先這樣建一個數值矩陣:

        x = -5:5
        t = (x.*x)' * ones(1,11)

然後這樣畫出來:

        [octave 語法]                   [rlab 語法]
        gsplot(t-t')                    splot(t-t')

這裡順便大力推薦 gnuplot, 實在超級好用, 理工科系的同學早學早賺, 不學會遺憾! 它本身的功能就只是畫圖而已, 既不能處理矩陣, 也不能做統計運算; 但這正是許多自由軟體共同的特點: 每套軟體只把自己的事情做好, 並提供與其他軟體良好的溝通管道 (例如其他軟體只要可以將數值整齊地印出來, gnuplot 就可以畫出來), 這樣學習者才能以簡馭繁, 發揮組合的力量。 [9]

我該學那一套?

筆者個人偏好 rlab, 因為它語言設計的「正交性」 (orthongonality) 比較好, 例如想取得矩陣 A 的行數, 在 rlab 裡面可以用 size(A)[2] (其實用 A.nc 更簡單) 但在 octave 裡面卻無法寫 size(A)(2) (不知道為什麼, 在命令列下可以; 但程式裡面就不行, 如果你幫我的 lademo.m 改成功了, 請告訴我)。 又如函式庫與個別函數是獨立的觀念, 不必有一對一的對應關係, 名稱也不必相同, 使用起來比較方便。 對於想寫 「作用於函數的 operators」 的 power user 來講, 函數可以當做參數來傳遞, 更是不可或缺的功能。 網站/手冊上開宗明義就說道 rlab 無意成為 matlab 的翻版, 而是要以更好的語法設計來取代它, 我相信 rlab 作者 Ian Searle 做到了。

另一方面, 用 octave 的人似乎比較多, 再加上 redhat 收錄的是 octave 而不是 rlab, 所以 octave 的聲勢似乎比較大。 由於與 matlab 相容, 因此既有的 matlab 工具箱也自然變成了 octave 使用者的資源。

不過自由軟體的世界與版權私有軟體不同, 我們認為多一些選擇總是好事, 就像 linux 有那麼多不同的版本, 更還有許多諸如 *BSD, Hurd 等等不同於 linux 的自由作業系統一樣。 畢竟這是一個尊重多元的時代, 不需要任何事情都定於一尊, 不是嗎?

參考網址

1. http://www.octave.org/
2. http://rlab.sourceforge.net/ 及 http://rlabplus.sourceforge.net/
3. http://www.r-project.org/
4. http://www-fourier.ujf-grenoble.fr/~parisse/giac.html
5. http://www.cyut.edu.tw/~ckhung/b/re/
6. http://www.octave.org/FAQ.html#SEC24
7. http://matlinks.sourceforge.net/
8. http://www.cyut.edu.tw/~ckhung/b/ma/gnuplot.php
9. http://www.cyut.edu.tw/~ckhung/b/018learn.php#combination
10. http://libai.math.ncu.edu.tw/bcc16/B/matlab/index.php 單維彰教授寫的 matlab 講義, 應該有一大半適用於 octave。

附錄: 正文當中的範例, 用 R 來打

R 的主要功能其實是統計運算. 如果你想從那個角度入門, 請參考 英文文件。 這裡我們只介紹它的向量與矩陣功能:

        # [快速導遊]
        3 * (7 - 2)
        x <- 3 * (7 - 2)     # 一旦指定入變數, 就不再顯示在螢幕上
        1/x
        x <- 3*(7-2); y <- 6.2e-3; z <- 4.3/8;
        sin(pi/3)
        sqrt(3)/2
        log(9)
        sqrt(3+4i)
        1/0
        0/0
        ? sqrt                  # 求助
        q()                     # 離開. 用 Ctrl-D 也可以.
        
        # [使用向量與矩陣]
        # 在 R 裡面, 向量是矩陣的基礎: 要建立矩陣之前, 必須先建立向量.
        r1 <- c(1,3,-13,4)   # 先建立三個列向量
        r2 <- c(5,-6,0,1)
        r3 <- c(-1,-7,5,2)
        a  <- rbind(r1,r2,r3) # 再將它們組成一個陣列
        # 或是建立一個很長的向量, 然後用 matrix 函數將它轉換成一個 3x4 的矩陣
        # 注意: column major 的順序!
        a  <- matrix(c(1,5,-1,3,-6,-7,-13,0,5,4,1,2), 3, 4)
        b  <- c(-3,5,-19,7)
        c  <- t(b)           # 轉置矩陣
        b; c                    # 注意: R 的向量, 內定視為行向量
        t(a)                    # 看看 「轉置」 作用在矩陣上的效果, 比較清楚.
        a[3,2]
        b[3]
        c[4]
        a[7]
        x  <- matrix(runif(48,0,1), 8, 6)
        x  <- floor(x*10)
        x[5,]
        x[,2]
        x[5:8,1:3]
        5:8
        seq(50,80,by=5)
        x[c(12,1:4),c(6,5)]
        matrix(c(x,x), 8, 12)
        t(matrix(c(t(x),t(x)), 6, 16))
        size(x)
        size(x)[1]
        size(x)[2]
        as.vector(a)
        t(as.vector(t(a)))
        sum(a)                  # R 的 sum 作用於整個陣列
        apply(a,2,sum)          # 若要分別對各行作用, 必須呼叫 apply
        apply(a,1,sum)          # 或想分別對各列作用也是.
        
        # [程式實例: Gauss-Jordan 消去法求 reduced row-echelon form]
        ...

---

• 本頁最新版網址: https://frdm.cyut.edu.tw/~ckhung/b/ma/linalg.php; 您所看到的版本: February 14 2012 02:32:25.
• 作者: 朝陽科技大學 資訊管理系 洪朝貴
• 寶貝你我的地球, 請 減少列印, 多用背面, 丟棄時做垃圾分類。
• 本文件以 Creative Commons Attribution-ShareAlike License 或以 Free Document License 方式公開授權大眾自由複製/修改/散佈。