该文介绍了java colt和commons-math3的一些矩阵计算API,并且使用colt库简单实现了基于法方程组法的最小二乘法,结构方程模型的梯度下降参数估计,广义混合效应模型(多层广义线性模型)的MCMC参数估计,实现和测试代码链接inuyasha11/stats
java矩阵计算概况
因为项目迁移需求,需要用java编写一些统计计算库。上网搜索了几个java矩阵库,找到了两个主流的,colt和commons-math3,colt库是CERN(欧洲核子研究组织)主导开发的,上次更新好像是10年前(?),所幸代码能支持java8,commons-math3一看名字就知道,系出apache软件基金会,里面除了矩阵库,还有其他的数学和统计方法,例如kmeans,遗传算法。这个两个库尝试了之后发现,真是难用,除了因为java缺乏操作符重载外,这两库API太少,很多轮子需要自己造,同时单线程速度也比底层调用blas,lpack的那些C++、python库慢,唯一的优点是纯java,可以很方便的多线程开发和后端系统对接,但是现在都搞微服务,谁还把计算模块嵌入业务模块啊。不过既然尝试,咱也写写使用介绍和心得吧,重点介绍下colt,稍微介绍下commons-math3。
矩阵构造
colt和commons-math3均支持两种形式的矩阵构造,dense矩阵和稀疏矩阵
对于dense矩阵,colt的实现是一维数组,commons-math3是二维数组,当矩阵元素大于4096的时候,commons-math3的dense矩阵工厂方法会创建BlockRealMatrix实例,BlockRealMatrix顾名思义,就是将矩阵分块存储,所以BlockRealMatrix虽然也是二维数组实现,但是数组的第一个索引仅仅是分块矩阵索引,第二个索引才是矩阵元素索引,本质上变成了和colt一样的一维数组实现。
对于稀疏矩阵,colt提供了稀疏矩阵类SparseDoubleMatrix2D,其存储矩阵元素的elements属性是一个hashmap,为了节省内存,colt自己实现了一个基于开放寻址方法的hashmap,在这个hashmap中,键为int类型,是矩阵的元素索引,值为double类型,是矩阵的元素值。当对稀疏矩阵set元素值的时候,如果元素值为0,则从hashmap中删除该元素的索引,若不为0,则在hashmap中存储键值对(索引,元素值),commons-math3的稀疏矩阵实现形式与colt类似,就不多赘述了。
colt和commoms-math3均支持通过传入矩阵的行数和列数构造初始化元素为0的矩阵
colt代码
import
commons-math3代码
RealMatrix
也可以通过传入二维数组构造矩阵
colt代码
double
commons-math3代码
double
也可以通过工厂模式创建矩阵实例
colt代码
DoubleFactory2D
这个工厂类还包含一些很有趣的静态方法,例如创建随机元素矩阵的方法random,合并两个矩阵的appendColumns方法和appendRows方法
DoubleFactory2D
commons-math3代码
RealMatrix
commons-math3的工厂方法会自动根据矩阵元素的大小创建Array2DRowRealMatrix类实例或BlockRealMatrix类实例,关于这两个类前面已经说明了一些情况,不多赘述了。
赋值和索引矩阵元素
colt
colt提供了get方法和getQuick方法读取矩阵中的单个元素值,还提供了set和setQuick方法设置矩阵中的单个元素值,get方法和getQuick方法的区别是get方法里加入了一些参数检查代码,set方法和setQuick方法同理,理论上getQuick和setQuick更快。
double
commons-math3
commons-math3提供了getEntry和setEntry两个方法索引和赋值矩阵元素,内部实现和colt几乎一样
double
获得矩阵的行数和列数
colt代码
int
commons-math3代码
int
复制矩阵
colt
在colt中,复制矩阵,只复制矩阵的形状,不复制元素,可以使用like方法
double
又复制矩阵的形状,又复制元素,可以使用copy方法
double
commons-math3
commons-math3貌似只提供了copy方法,commons-math3的copy方法调用了System.arraycopy静态方法进行数组赋值
RealMatrix
矩阵赋值速度比较
运行100次取平均值,单位毫秒
矩阵乘法
colt
colt矩阵的乘法可以通过zMult这个方法实现,zMult方法第1个参数为相乘矩阵实例,第2个参数为矩阵相乘的结果
double
如果zMult的第2个参数为null值,则zMult方法会自动为我们创建并返回结果矩阵
double
commons-math3
commons-math3提供了multiply方法进行矩阵乘法,multiply方法内部创建了新的矩阵实例并返回,所以不需要像colt那样传入储存结果得矩阵参数
RealMatrix
矩阵乘法速度比较
运行10次取平均值,单位毫秒,colt矩阵乘法运行时间是commons-math3的2倍
矩阵和矩阵的逐元素计算
colt
矩阵和矩阵的逐元素计算(例如矩阵的相减,相加,逐元素相乘,逐元素相除等)可以通过对assign方法传入另一个矩阵实例和cern.jet.math.Functions类的静态方法或自己实现DoubleDoubleFunction接口的类实例
import
assign方法的本质是两个嵌套循环,所以使用assign方法还是很慢的。
值得注意的是assign方法并不会创建新的矩阵实例,而是会在原矩阵实例基础上直接修改元素值,并返回原实例,这样做的好处是节省内存,但是会对代码编写造成很大的困扰,因为一不留神你传入的参数值就因为assign方法改变了,并且final关键字对于矩阵元素的修改并不起作用。
自己实现DoubleDoubleFunction接口类实例(下面用lambda方法简写了)作为assign方法的入参,下面的例子常在类似MCMC计算中出现,实现的计算是
double
类似python代码如下
mat1
commons-math3
commons-math3没有提供像colt的assign方法,只实现了基础的矩阵相加相减,并且commons-math矩阵相加相减都创建新的矩阵实例,一定程度会影响运行速度
// 加法
实现类似colt的嵌套循环assign并不是什么难事,只是速度很慢,例如自己通过commons-math3的getEntry和seEntry实现的类似colt assign的add方法(如下),运行速度比官方add方法慢
public
运行速度比较
运行10次取平均值,单位毫秒
矩阵和标量运算
colt
矩阵和标量的运算也是通过assign方法,下面代码实现的计算是
double
自己实现DoubleFunction接口的类实例
double
上面两段代码实现效果一样
commons-math3
略
矩阵的转置
colt
矩阵的转置可以通过viewDice方法实现,并且会创建一个新的矩阵实例并返回该新实例
double
commons-math3
略
矩阵的切片
略
线性代数
colt
colt提供了一些基础的线性代数算法,cern.colt.matrix.linalg.Algebra实例提供了解方程方法,
CholeskyDecomposition实例提供了cholesky分解,还有svd分解,qr分解等,详细可以参考应用实战的最小二乘法
commons-math3
略
统计
colt
cern.colt.matrix.doublealgo.Statistic类提供了一些静态方法计算矩阵的方差协方差矩阵,相关系数
DoubleFactory2D
commons-math3
略
应用实战
应用实战使用colt做实例
最小二乘之法方程组法
法方程组发实现细节如下
实现目标,自变量,响应变量,特征
计算
计算
计算cholesky分解,
解方程组,
转换为代码部分如下
/**
法方程组法主要用到了colt提供的CholeskyDecomposition类进行cholesky分解和Algebra类进行解方程组
结构方程模型之一阶验证性因子分析
关于结构方程模型细节,可以参考作者专栏写的另外三篇文章
这里实现的是最简单的结构方程模型,一阶验证性因子分析,固定因子方差为1,设定误差之间相互独立
colt并没有提供多元正态分布随机抽样方法,这边为了测试代码,还需要制造模拟数据,而模拟数据需要用到多元正态分布随机抽样。所幸,colt提供了一元正态分布随机采样和SVD分解,通过一元正态分布随机采样和SVD分解我们可以随机采样多元正态分布,下面代码实现的是均值为0向量的多元正态分布随机采样,整体步骤是(1)随机采样标准正态分布(2)计算方差协方差的svd分解(3)一系列矩阵计算得到多元正态分布的随机数。
/**
参数估计
一阶验证性因子分析参数估计过程
实现目标,样本协方差矩阵,
约束对角线元素为1,除了对角线元素其余为0,初值为0的元素恒为0
输入观察变量,初值,步长,最大迭代次数
计算的有偏方差协方差矩阵
计算,,的梯度
依据梯度和步长更新,,
重复上述两个步骤直到达到最大迭代次数
梯度下降的部分代码如下
/**
广义混合效应模型(多层广义线性模型)之MCMC参数估计
这里的的广义混合效应模型是最简单的广义混合效应模型,关于混合效应模型,可以参考作者专栏写的另一篇文章
换一种写法
我们通过MCMC吉布斯采样进行参数估计,整体过程如下
(1) 输入响应变量,链长,burn,thin,设,
(2),的设为0向量
(3)随机采样,
(4)计算后验似然函数与后验似然函数的比值
(5)随机采样,
(6) 若,
(7)随机采样,
(8)计算后验似然函数与后验似然函数的比值
(9)随机采样,
(10) 若,
(11)重复(3)到(10)直至迭代结束
部分代码如下
首先,我们需要实现一个logistic函数
/**
注意到一个小技巧,为了防止数值溢出,作者给logistic值设定了1个上界和下界,这是数值分析常用的做法。
colt 并没有实现不同形状向量相加方法,导致我们不得不写一个类似于外积的外和方法,主要是求
/**
完整的实现和测试代码地址如下
inuyasha11/statsgitee.com