该文介绍了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/stats​gitee.com