Octave：矩阵计算的新宠

 

作者：于江生（北京大学计算机系）

声明：允许未经作者的同意进行非商业目的的转载，但必须保持原文的完整性。

---

实话实说，MatLab是迄今为止矩阵计算最强大的工具（没有之一）。可惜MatLab是商用的，一般个体还真买不起。MatLab的Windows版本比Linux版本要好些，这让我不敢轻易断言Windows一无是处，毕竟其下有MatLab这样强悍的软件。以前在Windows下工作，MatLab一直是我的首选矩阵计算工具，在统计计算工具S-PLUS出现之前，人们快乐地用着MatLab简陋的统计工具箱。后来有了R，它彻底地坐稳了统计计算的头把交椅，MatLab似乎也无意去争夺全料冠军，但事实上它在很多方面都做得无可挑剔。这让我们这些买不起却很需要MatLab的穷人感慨不已，MatLab如果是免费的该多好……

为何选择使用octave？

导入文件

Octave与MatLab的一些小区别

布尔值的乘积

逻辑运算符、算术运算符

C-风格的自动增量、赋值、屏幕打印

注意空格

直方图内置函数hist

导入空文件

行续符

if、for等环境的结束符

R和octave命令的对照表

R读入octave导出的数据

为何选择使用octave？

SciLab和octave是开源的且免费的矩阵计算工具，二者都有希望成为矩阵计算的新宠。相比之下，

• octave与MatLab的兼容性更高。
• octave遵循GPL协议（GNU General Public License），用户可以单独发行octave或者包含在其产品中发行。而scilab则不允许，你只能免费地使用它。
• octave没有图形界面，是命令交互的。在某些人眼里这是不可饶恕的缺点，而在另外一些人眼里则是大大的优点。

它们都具备以下特点：以矩阵为基本数据类型，内置支持复数，有内置函数和外部函数库，用户自定义函数的可扩展性等特点。UNIX的很多用户选择使用octave，看中的就是它与MatLab兼容性好这一事实。随着开源运动的深入人心，octave不断地发展壮大，它会吸引一大批MatLab的使用者。

GNU octave网站：http://www.octave.org/

好习惯从头开始：

• 首先学会使用help，搞不定再到网上查，最后才求人。
• 学习octave的捷径：读octave的函数源码。
• 每个命令都以“;”结束，否则矩阵的具体内容会显示出来。
• 学会适当地使用内置命令clear，从内存中清除一些无用数据或变元。
• 如果没有必要，不要轻易改变矩阵大小。
• 重要的中间结果要保存。

---

导入文件

octave和MatLab一样用load导入数据文件，譬如

octave>  A = load data.txt ;

将把data.txt里的数据导入octave并赋给矩阵 A。对于图像文件，octave用imread将图像导入并存为矩阵img，

img = imread("jam.jpg") ;

在octave里显示图像很简单，用命令:

imshow(img) ;

除了jpeg和png格式的图像可以直接导入，其他格式的图像必须经过ImageMagick的convert函数转换后才可读入。ImageMagick是命令行的强大的图像处理工具，convert几乎涵盖了所有格式图像的转换。

如果你关心imread函数的源码，可以去读 /usr/local/share/octave/packages/image-1.0.8/imread.m，该函数把灰度图像导入为MxN矩阵，把彩色图像导入为MxNx3矩阵。具体的帮助文件，可以

help imread ;

或者来个更详细点儿的

help -i imread ;

---

Octave与MatLab的一些小区别

MatLab用户转而使用octave几乎不需要什么培训，只是要一些小细节上注意一下。下面我们罗列一些octave和MatLab的区别。

布尔值的乘积

X = ones(2,2) ;prod(size(X)==1)

MatLab和octave的输出是不同的：

Matlab: ??? Function 'prod' is not defined for values of class 'logical'.Octave: ans = 0

octave输出为0的原因是 size(X) 为

ans =   2   2

逻辑运算符、算术运算符

Octave与MatLab兼容，甚至更为宽松。如，

运算Matlaboctave或|“|” 或者“||”且&& 或者 &&否~=~= 或者 !=

MatLab用 x^2，octave用 x^2 或者 x**2 表示 “x的平方”。Octave用 x**2 是为了照顾GnuPlot的用户。总而言之，octave在运算符方面彻底兼容MatLab，MatLab用户放心大胆地用octave吧，但octave用户用MatLab的时候就要小心了。

C-风格的自动增量、赋值、屏幕打印

Octave允许C-风格的

i++ ; ++i ; i+=1 ;printf('My result is: %d/n', 4)

而MatLab不认它们。MatLab打印至屏幕和文件都用 fprintf 函数。

注意空格

octave对空格是作为一个符号识别的，在列合并中短的列自然扩充，例如

A = ['123 ';'123'] ;size(A)

的结果是 2 4，而MatLab则返回列合并有问题：

?? Error using ==> vertcat

另外，转置符号与矩阵之间如果有空格

[0 1] '

在MatLab里不允许，octave则允许，且与 [0 1]' 的结果是一样的。

直方图内置函数hist

octave的hist为

hist (Y, X, NORM)

其中NORM为所有柱高之和。

导入空文件

MatLab允许导入空文件，老版本的octave不允许，新版本的octave-3.0.3则允许。

行续符

MatLab中用 `...' 做行续符，如用

A = rand (1, ...         2) ;

表达

A = rand (1,2) ;

Octave与MatLab兼容，除此之外，octave还允许如下两种表示方法。

A = rand (1,        2) ;

和

A = rand (1, /        2) ;

if、for等环境的结束符

Octave用

end{if,for, ...}

而MatLab则统一用 end。

---

R和octave命令的对照表

octave和R联合起来用的时候，我们需要下面的命令对照表帮助我们理清楚它们的区别。“无”仅仅是说没有一个命令行的简单表示，并不代表不能表示。这种对比不是比较谁更强大，而是为了记忆，无论是对R用户学习octave或者octave用户学习R，都是有所裨益的。

octaveR帮助help -ihelp.start()helphelp(help)help sorthelp(sort)demo()lookfor plotapropros('plot')help.search('plot')复数3+4i3+4ii1i % R把"i"视为变量名abs(3+4i)Mod(3+4i)arg(3+4i)Arg(3+4i)conj(3+4i)Conj(3+4i)real(3+4i)Re(3+4i)imag(3+4i)Im(3+4i)向量、序列1:101:10 或 seq(10)1:3:10seq(1,10,by=3)10:-1:110:110:-3:1seq(from=10,to=1,by= -3)linspace(1,10,7)seq(1,10,length=7)(1:10)+i1:10+1ia=[2 3 4 5]; # 不显示结果a <- c(2,3,4,5) % 不用加分号a=[2 3 4 5] #显示结果(a <- c(2,3,4,5)) % 显示结果adash=[2 3 4 5]' ;adash <- t(c(2,3,4,5))[a a]c(a,a)[a a*3]c(a,a*3)a.*aa*aa.^3a^3向量的合并与重复[1:4 a]c(1:4,a)[1:4 1:4]rep(1:4,2)无rep(1:4,1:4) % 结果是：1 2 2 3 3 3 4 4 4 4无rep(1:4,each=3) % 结果是：1 1 1 2 2 2 3 3 3 4 4 4a=1:100;a <- 1:100a(2:100)a[-1] % a 去掉第1个元素a([1:9 11:100])a[-10] % a 去掉第10个元素无a[-seq(1,50,3)] % a 去掉第1,4,7,...个元素向量的赋值a(a>90)= -44;a[a>90] <- -44向量的最大、最小a=randn(1,4);a <- rnorm(4)b=randn(1,4);b <- rnorm(4)max(a,b)pmax(a,b)max([a' b'])cbind(max(a),max(b))max([a b])max(a,b)[m i] = max(a)m <- max(a) ; i <- which.max(a)"min" 类似向量的秩ranks(rnorm(8,1))rank(rnorm(8))ranks(rnorm(randn(5,6)))apply(matrix(rnorm(30),6),2,rank)矩阵的行合并与列合并[1:4 ; 1:4]rbind(1:4,1:4)[1:4 ; 1:4]'cbind(1:4,1:4) 或 t(rbind(1:4,1:4))[2 3 4 5]c(2,3,4,5)[2 3;4 5]rbind(c(2,3),c(4,5)) % rbind() 合并行; cbind() 合并列[2 3;4 5]'cbind(c(2,3),c(4,5)) 或 matrix(2:5,2,2)a=[5 6];a <- c(5,6)b=[a a;a a];b <- rbind(c(a,a),c(a,a))[1:3 1:3 1:3 ; 1:9]rbind(1:3, 1:9)[1:3 1:3 1:3 ; 1:9]'cbind(1:3, 1:9)无rbind(1:3, 1:8)产生矩阵ones(4,7)matrix(1,4,7) 或 array(1,c(4,7))ones(4,7)*9matrix(9,4,7) 或 array(9,c(4,7))eye(3)diag(1,3) % 对角线都为1的对角阵diag([4 5 6])diag(c(4,5,6)) % 对角线为4,5,6的对角阵diag(1:10,3)无reshape(1:6,2,3)matrix(1:6,nrow=2) 或 array(1:6,c(2,3))reshape(1:6,3,2)matrix(1:6,ncol=2) 或 array(1:6,c(3,2))reshape(1:6,3,2)'matrix(1:6,nrow=2,byrow=T)a=reshape(1:36,6,6);a <- matrix(1:36,c(6,6))rem(a,5)a %% 5a(rem(a,5)==1)= -999a[a%%5==1] <- -999a(:)as.vector(a)矩阵中抽取元素a=reshape(1:12,3,4);a <- matrix(1:12,nrow=3)a(2,3)a[2,3]a(2,:)a[2, ]a(2:3,:)a[-1,]a(:,[1 3 4])a[,-2]a(:,1)a[ ,1]a(:,2:4)a[ ,-1]a([1 3],[1 2 4])a[-2,-3]矩阵赋值a(:,1) = 99a[ ,1] <- 99a(:,1) = [99 98 97]'a[ ,1] <- c(99,98,97)矩阵：转置、共轭a'Conj(t(a))a.'t(a)矩阵：求和a=ones(6,7)a <- matrix(1,6,7)sum(a)apply(a,2,sum)sum(a')apply(a,1,sum)sum(sum(a))sum(a)cumsum(a)apply(a,2,cumsum)cumsum(a')apply(a,1,cumsum)矩阵排序a=rand(3,4);a <- matrix(runif(12),c(3,4))sort(a(:))sort(a)sort(a)apply(a,2,sort)sort(a')apply(a,1,sort)cummax(a)apply(a,2,cummax)矩阵：最大、最小a=randn(100,4)a <- matrix(rnorm(400),4)max(a)apply(a,1,max)[v i] = max(a)v <- apply(a,1,max) ; i <- apply(a,1,which.max)b=randn(4,4);b <-matrix(rnorm(16),4)c=randn(4,4);c <-matrix(rnorm(16),4)max(b,c)pmax(b,c)矩阵的乘法a=reshape(1:6,2,3);a <- matrix(1:6,2,3)b=reshape(1:6,3,2);b <- matrix(1:6,3,2)c=reshape(1:4,2,2);c <- matrix(1:4,2,2)v=[10 11];v <- c(10,11)w=[100 101 102];w <- c(100,101,102)x=[4 5]' ;x <- t(c(4,5))a*ba %*% bv*av %*% aa*w'a %*% wb*v'b %*% vv*xx %*% v 或 v %*% t(x)x*vt(x) %*% vv*a*w'v %*% a %*% wv .* x'v*x 或_ x*va .* [w ;w]w * aa .* [x x x]a * t(rbind(x,x,x)) 或 a*as.vector(x)v*cv %*% cc*v'c %*% v其他矩阵操作a=rand(3,4);a <- matrix(runif(12),c(3,4))fliplr(a)a[,4:1]flipud(a)a[3:1,]a=reshape(1:9,3,3)a <- matrix(1:9,3)vec(a)as.vector(a)vech(a)a[row(a) <= col(a)]size(a)dim(a)网格[x y]=meshgrid(1:5,10:12);无查找find(1:10 > 5.5)which(1:10 > 5.5)a=diag([4 5 6])a <- diag(c(4,5,6))find(a)which(a != 0) % which() 的变元是布尔变元[i j]= find(a)which(a != 0,arr.ind=T)[i j k]=find(a)ij <- which(a != 0,arr.ind=T); k <- a[ij]读文件load foo.txtf <- read.table("~/foo.txt")f <- as.matrix(f)写文件save -ascii bar.txt fwrite(f,file="bar.txt")图形输出gset output "foo.eps"postscript(file="foo.eps")gset terminal postscript epsplot(1:10)plot(1:10)dev.off ()赋值string="a=234";string <- "a <- 234"eval(string)eval(parse(text=string))产生随机数均匀分布rand(10,1)runif(10)2+5*rand(10,1)runif(10,min=2,max=7) 或 runif(10,2,7)rand(10)matrix(runif(100),10)正态分布randn(10,1)rnorm(10)2+5*randn(10,1)rnorm(10,2,5)rand(10)matrix(rnorm(100),10)beta分布hist(beta_rnd(4,2,1000,1)hist(rbeta(1000,shape1=4,shape2=10)) 或 hist(rbeta(1000,4,10))FOR循环for i=1:5; disp(i); endforfor(i in 1:5) {print(i)}多项式的根roots([1 2 1])polyroot(c(1,2,1))polyval([1 2 1 2],1:10)无集合论a = create_set([1 2 2 99 2 ])a <- sort(unique(c(1,2,2,99,2)))b = create_set([2 3 4 ])b <- sort(unique(c(2,3,4)))intersect(a,b)intersect(a,b)union(a,b)union(a,b)complement(a,b)setdiff(b,a)any(a == 2)is.element(2,a)绘图a=rand(10);a <- array(runif(100),c(10,10))help plothelp (plot) and methods(plot)plot(a)matplot(a,type="l",lty=1)plot(a,'r')matplot(a,type="l",lty=1,col="red")plot(a,'x')matplot(a,pch=4)plot(a,'—')matplot(a,type="l",lty=2)plot(a,'x-')matplot(a,pch=4,type="b",lty=1)plot(a,'x—')matplot(a,pch=4,type="b",lty=2)semilogy(a)matplot(a,type="l",lty=1,log="y")semilogx(a)matplot(a,type="l",lty=1,log="x")loglog(a)matplot(a,type="l",lty=1,log="xy")plot(1:10,'r')plot(1:10,col="red",type="l")hold onmatplot(10:1,col="blue",type="l",add=T)plot(10:-1:1,'b')gridgrid()a=randn(10);a <- matrix(rnorm(100),nr=10)contour(a)contour(a)contour(a,77)contour(a,nlevels=77) ; filled.contour(a)mesh(rand(10))persp(matrix(runif(100),10),theta=30,phi=30,d=1e9)文件与操作系统system("ls")system("ls")pwdgetwd()cdsetwd()

R读入octave导出的数据

统计计算软件R的 foreign 包提供了函数 read.octave，可以读入 octave 用命令 save -ascii 创建的文本数据文件，且支持变量的大多数通用类型，包括标准的原子型（复矩阵， N维数组，字符串，布尔矩阵等）和 递归式（结构体，单元和列表）。

在octave中用 save -ascii 保存的矩阵数据，也可以在 R 中用命令 read.table 导入，然后用 as.matrix() 强制为 R 中的矩阵使用。我比较倾向于这种方法。这样我们就能充分利用octave擅长矩阵计算和R擅长统计计算的优势，将二者联合起来使用。我们将详细介绍octave读入图像文件，输出能被R处理的矩阵数据。

下面举个例子：我们投掷一枚硬币，已知正面出现的概率为 p，恰好掷出 R 正面所用的次数 N 是我们要考察的，我们做 E 次随机试验，看看N的经验分布情况。

## File      : toss.m## Purpose   : The numbers of tossing to get R heads## Author    : Jiangsheng Yu ([e-mail])## Data      : 11-26-2008## Available : http://icl.pku.edu.cn/member/yujs/Computing.htm## Usage     : run toss.mmore off ;   ## turn the pagination offE = 10000;   ## the number of experimentsresult = zeros(E,1);   ## the sequence of E resultsR = 6   ;    ## required number of headsp = 0.3 ;    ## the probability of headH = 0 ;      ## no heads at the beginningN = 0 ;      ## no tosses at the beginningfor i = 1:E    do  ## if head, outcome=1; otherwise, outcome=0  outcome = (rand(1,1) < p) ;  H += outcome ; ## the total number of heads  N += 1 ;       ## the total number of tossesuntil ( H >= R ) ## until R headsresult(i,1) = N ;N = 0 ;H = 0 ;endforhist (result,40,1) ;

对于 p=0.3 ，R=2 做 E=10000 次随机试验得到 N 的直方图如下：掷出R个正面所用次数的直方图

对于 p=0.3 ，R=6 做 E=10000 次随机试验得到 N 的直方图如下：掷出R个正面所用次数的直方图

我们把结果保存为 result.data，再读到 R 中处理这些数据。

> x = result' ;> save result.dat x ;

在 R 中我们读入数据，然后画出直方图。

> library(foreign)> a <- read.octave("result.dat")> hist(a$x, freq= FALSE, col="blue", border="pink")

得到 p=0.3 ，R=6 的直方图：掷出R个正面所用次数的直方图