牛厨馆 / 待分类 / Octave:矩阵计算的新宠

分享

   

Octave:矩阵计算的新宠

2020-07-04  牛厨馆
                                                

作者:于江生(北京大学计算机系)

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


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

为何选择使用octave?

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

  • octave与MatLab的兼容性更高。

  • octave遵循GPL协议(GNU General Public License),用户可以单独发行octave或者包含在其产品中发行。而scilab则不允许,你只能免费地使用它。

  • octave没有图形界面,是命令交互的。在某些人眼里这是不可饶恕的缺点,而在另外一些人眼里则是大大的优点。

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

GNU octave网站:

好习惯从头开始:

  • 首先学会使用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+4i
i1i % 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:1
10:-3:1seq(from=10,to=1,by= -3)
linspace(1,10,7)seq(1,10,length=7)
(1:10)+i1:10+1i
a=[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*a
a.^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 4
a=1:100;a <- 1:100
a(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 %% 5
a(rem(a,5)==1)= -999a[a%%5==1] <- -999
a(:)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] <- 99
a(:,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 %*% b
v*av %*% a
a*w'a %*% w
b*v'b %*% v
v*xx %*% v 或 v %*% t(x)
x*vt(x) %*% v
v*a*w'v %*% a %*% w
v .* x'v*x 或_ x*v
a .* [w ;w]w * a
a .* [x x x]a * t(rbind(x,x,x)) 或 a*as.vector(x)
v*cv %*% c
c*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 (yujs@pku.edu.cn)
## Data      : 11-26-2008
## Available : http://icl.pku.edu.cn/member/yujs/Computing.htm
## Usage     : run toss.m

more off ;   ## turn the pagination off
E = 10000;   ## the number of experiments
result = zeros(E,1);   ## the sequence of E results
R = 6   ;    ## required number of heads
p = 0.3 ;    ## the probability of head
H = 0 ;      ## no heads at the beginning
N = 0 ;      ## no tosses at the beginning
for 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 tosses
	until ( H >= R ) ## until R heads
	result(i,1) = N ;
	N = 0 ;
	H = 0 ;
endfor
hist (result,40,1) ;

对于 p=0.3 ,R=2 做 E=10000 次随机试验得到 N 的直方图如下:

掷出R个正面所用次数的直方图
掷出R个正面所用次数的直方图

对于 p=0.3 ,R=6 做 E=10000 次随机试验得到 N 的直方图如下:

掷出R个正面所用次数的直方图
掷出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个正面所用次数的直方图
掷出R个正面所用次数的直方图

    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多
    喜欢该文的人也喜欢 更多

    ×
    ×

    ¥.00

    微信或支付宝扫码支付:

    开通即同意《个图VIP服务协议》

    全部>>