该文介绍了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 cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
DoubleMatrix2D mat = new DenseDoubleMatrix2D(3, 4);
System.out.println(mat.toString());
commons-math3代码
RealMatrix mat = new Array2DRowRealMatrix(4,4);
也可以通过传入二维数组构造矩阵
colt代码
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}, {7d, 8d, 9d}};
DoubleMatrix2D mat = new DenseDoubleMatrix2D(matData);
System.out.println(mat.toString());
commons-math3代码
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}, {7d, 8d, 9d}};
RealMatrix mat = new Array2DRowRealMatrix(matData);
也可以通过工厂模式创建矩阵实例
colt代码
DoubleFactory2D F = DoubleFactory2D.dense;
F.make(4,4);
这个工厂类还包含一些很有趣的静态方法，例如创建随机元素矩阵的方法random，合并两个矩阵的appendColumns方法和appendRows方法
DoubleFactory2D F = DoubleFactory2D.dense;
DoubleMatrix2D mat1 = F.random(10, 5);
DoubleMatrix2D mat2 = F.random(10, 3);
DoubleMatrix2D mat3 = F.appendColumns(mat1, mat2);
System.out.println(mat3);
commons-math3代码
RealMatrix mat = MatrixUtils.createRealMatrix(4, 4);
commons-math3的工厂方法会自动根据矩阵元素的大小创建Array2DRowRealMatrix类实例或BlockRealMatrix类实例，关于这两个类前面已经说明了一些情况，不多赘述了。
赋值和索引矩阵元素
colt
colt提供了get方法和getQuick方法读取矩阵中的单个元素值，还提供了set和setQuick方法设置矩阵中的单个元素值，get方法和getQuick方法的区别是get方法里加入了一些参数检查代码，set方法和setQuick方法同理，理论上getQuick和setQuick更快。
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}, {7d, 8d, 9d}};
DoubleMatrix2D mat = new DenseDoubleMatrix2D(matData);
System.out.println(mat.get(0, 0));
mat.set(0, 0, 10d) ;
System.out.println(mat.get(0, 0));
commons-math3
commons-math3提供了getEntry和setEntry两个方法索引和赋值矩阵元素，内部实现和colt几乎一样
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}, {7d, 8d, 9d}};
RealMatrix mat = new Array2DRowRealMatrix(matData);
double val = mat.getEntry(0, 0);
mat.setEntry(0, 0, val + 10d);
获得矩阵的行数和列数
colt代码
int rowSize = mat.rows();
int colSize = mat.rows();
commons-math3代码
int rowSize = mat.getRowDimension();
int colSize = mat.getColumnDimension();
复制矩阵
colt
在colt中，复制矩阵，只复制矩阵的形状，不复制元素，可以使用like方法
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}};
DoubleMatrix2D mat1 = new DenseDoubleMatrix2D(matData);
DoubleMatrix2D mat2 = mat1.like();
System.out.println(mat2);
又复制矩阵的形状，又复制元素，可以使用copy方法
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}};
DoubleMatrix2D mat1 = new DenseDoubleMatrix2D(matData);
DoubleMatrix2D mat2 = mat1.copy();
System.out.println(mat2);
commons-math3
commons-math3貌似只提供了copy方法，commons-math3的copy方法调用了System.arraycopy静态方法进行数组赋值
RealMatrix mat2 = mat1.copy();
矩阵赋值速度比较
运行100次取平均值，单位毫秒
矩阵乘法
colt
colt矩阵的乘法可以通过zMult这个方法实现，zMult方法第1个参数为相乘矩阵实例，第2个参数为矩阵相乘的结果
double[][] matData1 = {{1d, 2d, 3d}, {4d, 5d, 6d}};
double[][] matData2 = {{1d, 2d}, {3d, 4d}, {5d, 6d}};
DoubleMatrix2D mat1 = new DenseDoubleMatrix2D(matData1);
DoubleMatrix2D mat2 = new DenseDoubleMatrix2D(matData2);
// mat1乘mat2后得到的矩阵
DoubleMatrix2D mat3 = new DenseDoubleMatrix2D(mat1.rows(), mat2.columns());
mat1.zMult(mat2, mat3);
System.out.println(mat3);
如果zMult的第2个参数为null值，则zMult方法会自动为我们创建并返回结果矩阵
double[][] matData1 = {{1d, 2d, 3d}, {4d, 5d, 6d}};
double[][] matData2 = {{1d, 2d}, {3d, 4d}, {5d, 6d}};
DoubleMatrix2D mat1 = new DenseDoubleMatrix2D(matData1);
DoubleMatrix2D mat2 = new DenseDoubleMatrix2D(matData2);
// mat1乘mat2后得到的矩阵
DoubleMatrix2D mat3 = mat1.zMult(mat2, null);
System.out.println(mat3);
commons-math3
commons-math3提供了multiply方法进行矩阵乘法，multiply方法内部创建了新的矩阵实例并返回，所以不需要像colt那样传入储存结果得矩阵参数
RealMatrix mat3 = mat1.multiply(mat2);
矩阵乘法速度比较
运行10次取平均值，单位毫秒，colt矩阵乘法运行时间是commons-math3的2倍
矩阵和矩阵的逐元素计算
colt
矩阵和矩阵的逐元素计算（例如矩阵的相减，相加，逐元素相乘，逐元素相除等）可以通过对assign方法传入另一个矩阵实例和cern.jet.math.Functions类的静态方法或自己实现DoubleDoubleFunction接口的类实例
import cern.jet.math.Functions;
double[][] matData1 = {{1d, 2d, 3d}, {4d, 5d, 6d}};
double[][] matData2 = {{7d, 8d, 9d}, {10d, 11d, 12d}};
DoubleMatrix2D mat1 = new DenseDoubleMatrix2D(matData1);
DoubleMatrix2D mat2 = new DenseDoubleMatrix2D(matData2);
// mat1加mat2
mat1.assign(mat2, Functions.plus);
// mat1减mat2
mat1.assign(mat2, Functions.minus);
assign方法的本质是两个嵌套循环，所以使用assign方法还是很慢的。
值得注意的是assign方法并不会创建新的矩阵实例，而是会在原矩阵实例基础上直接修改元素值，并返回原实例，这样做的好处是节省内存，但是会对代码编写造成很大的困扰，因为一不留神你传入的参数值就因为assign方法改变了，并且final关键字对于矩阵元素的修改并不起作用。
自己实现DoubleDoubleFunction接口类实例（下面用lambda方法简写了）作为assign方法的入参，下面的例子常在类似MCMC计算中出现，实现的计算是
double[][] matData1 = {{1d, 2d, 3d}, {4d, 5d, 6d}};
double[][] matData2 = {{2d, 1d, 4d}, {5d, 3d, 2d}};
DoubleMatrix2D mat1 = new DenseDoubleMatrix2D(matData1);
DoubleMatrix2D mat2 = new DenseDoubleMatrix2D(matData2);
mat1.assign(mat2, (x, y) -> x < y ? x:y);
System.out.println(mat1);
类似python代码如下
mat1 = numpy.array([[1, 2, 3], [4, 5, 6]])
mat2 = numpy.array([[2, 1, 4], [5, 3, 2]])
mat1[mat1 > mat2] = mat2[mat1 > mat2]
commons-math3
commons-math3没有提供像colt的assign方法，只实现了基础的矩阵相加相减，并且commons-math矩阵相加相减都创建新的矩阵实例，一定程度会影响运行速度
// 加法
RealMatrix mat3 = mat1.add(mat2);
// 减法
RealMatrix mat4 = mat1.subtract(mat2);
实现类似colt的嵌套循环assign并不是什么难事，只是速度很慢，例如自己通过commons-math3的getEntry和seEntry实现的类似colt assign的add方法（如下），运行速度比官方add方法慢
public static void add(RealMatrix x, RealMatrix y) {
int rowSize = x.getRowDimension();
int colSize = x.getColumnDimension();
for (int i = 0; i < rowSize; i++) {
for (int j = 0; j < colSize; j++) {
x.setEntry(i, j, x.getEntry(i, j) + y.getEntry(i, j));
}
}
}
运行速度比较
运行10次取平均值，单位毫秒
矩阵和标量运算
colt
矩阵和标量的运算也是通过assign方法，下面代码实现的计算是
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}};
DoubleMatrix2D mat = new DenseDoubleMatrix2D(matData);
// 矩阵所有元素加0.5
mat.assign(Functions.plus(0.5));
System.out.println(mat);
自己实现DoubleFunction接口的类实例
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}};
DoubleMatrix2D mat = new DenseDoubleMatrix2D(matData);
// 矩阵所有元素加0.5
mat.assign(x -> x + 0.5);
System.out.println(mat);
上面两段代码实现效果一样
commons-math3
略
矩阵的转置
colt
矩阵的转置可以通过viewDice方法实现，并且会创建一个新的矩阵实例并返回该新实例
double[][] matData = {{1d, 2d, 3d}, {4d, 5d, 6d}};
DoubleMatrix2D mat = new DenseDoubleMatrix2D(matData);
DoubleMatrix2D matT = mat.viewDice();
System.out.println(matT.toString());
commons-math3
略
矩阵的切片
略
线性代数
colt
colt提供了一些基础的线性代数算法，cern.colt.matrix.linalg.Algebra实例提供了解方程方法，
CholeskyDecomposition实例提供了cholesky分解，还有svd分解，qr分解等，详细可以参考应用实战的最小二乘法
commons-math3
略
统计
colt
cern.colt.matrix.doublealgo.Statistic类提供了一些静态方法计算矩阵的方差协方差矩阵，相关系数
DoubleFactory2D F = DoubleFactory2D.dense;
DoubleMatrix2D mat = F.random(100, 5);
// 方差协方差矩阵
DoubleMatrix2D cov = Statistic.covariance(mat);
System.out.println(cov.toString());
commons-math3
略
应用实战
应用实战使用colt做实例
最小二乘之法方程组法
法方程组发实现细节如下
，
自变量，
响应变量，
特征
计算
计算
计算cholesky分解，
解方程组，
转换为代码部分如下
/**
* 自变量
*/
private final DoubleMatrix2D a;
/**
* 响应变量
*/
private final DoubleMatrix2D b;
/**
* 特征
*/
private DoubleMatrix2D x;
/**
* 法方程组法最小二乘法
* 参数估计
*/
public void solve() {
DoubleMatrix2D at = a.viewDice();
DoubleMatrix2D aa = at.zMult(a, null);
// cholesky分解
CholeskyDecomposition cholesky = new CholeskyDecomposition(aa);
DoubleMatrix2D l = cholesky.getL();
DoubleMatrix2D ab = at.zMult(b, null);
Algebra algebra = new Algebra();
DoubleMatrix2D temp = algebra.solve(l, ab);
x = algebra.solve(l.viewDice(), temp);
法方程组法主要用到了colt提供的CholeskyDecomposition类进行cholesky分解和Algebra类进行解方程组
结构方程模型之一阶验证性因子分析
关于结构方程模型细节，可以参考作者专栏写的另外三篇文章
这里实现的是最简单的结构方程模型，一阶验证性因子分析，固定因子方差为1，设定误差之间相互独立
colt并没有提供多元正态分布随机抽样方法，这边为了测试代码，还需要制造模拟数据，而模拟数据需要用到多元正态分布随机抽样。所幸，colt提供了一元正态分布随机采样和SVD分解，通过一元正态分布随机采样和SVD分解我们可以随机采样多元正态分布，下面代码实现的是均值为0向量的多元正态分布随机采样，整体步骤是（1）随机采样标准正态分布（2）计算方差协方差的svd分解（3）一系列矩阵计算得到多元正态分布的随机数。
/**
* 标准正态分布的随机矩阵采样静态方法
* @param rowSize 矩阵的行数
* @param colSize 矩阵的列数
* @return 标准正态分布随机矩阵
*/
public static double[][] genRandomStdNormData(int rowSize, int colSize) {
double[][] data = new double[rowSize][colSize];
for (int row = rowSize; --row >= 0; ) {
for (int col = colSize; --col >= 0; ) {
data[row][col] = Normal.staticNextDouble(0, 1);
}
}
return data;
}
/**
* 均值为0向量的多元正态分布随机采样
* @param size 样本量
* @param cov 方差协方差矩阵
* @return 多元正态分布随机采样样本
*/
public static DoubleMatrix2D genRandomMeanZeroMvnFactor(int size, DoubleMatrix2D cov) {
double[][] xData = genRandomStdNormData(size, cov.rows());
DoubleMatrix2D x = new DenseDoubleMatrix2D(xData);
// svd分解
SingularValueDecomposition svd = new SingularValueDecomposition(cov);
DoubleMatrix2D s = svd.getS();
DoubleMatrix2D v = svd.getV();
Algebra algebra = new Algebra();
DoubleMatrix2D vh = algebra.inverse(v);
DoubleMatrix2D sv = s.assign(Functions.sqrt).zMult(vh, null);
return x.zMult(sv, null);
}
参数估计
一阶验证性因子分析参数估计过程
，
样本协方差矩阵，
约束对角线元素为1，
除了对角线元素其余为0，
初值为0的元素恒为0
输入观察变量，
初值，步长，最大迭代次数
计算的有偏方差协方差矩阵
计算，
，
的梯度
依据梯度和步长更新，
，
重复上述两个步骤直到达到最大迭代次数
梯度下降的部分代码如下
/**
* 因子载荷
*/
private DoubleMatrix2D lambda;
/**
* 因子方差协方差矩阵
*/
private DoubleMatrix2D phi;
/**
* 误差方差协方差矩阵
*/
private DoubleMatrix2D error;
/**
* 样本有偏方差协方差矩阵
*/
private DoubleMatrix2D s;
/**
* 梯度下降步长
*/
private final double stepSize;
/**
* 梯度下降迭代次数
*/
private final int maxIter;
/**
* 参数梯度下降
*/
private void gradientDescent() {
DoubleMatrix2D sigma = lambda.zMult(phi, null).zMult(lambda.viewDice(), null);
sigma.assign(error, plus);
DoubleMatrix2D sigmaInv = new Algebra().inverse(sigma);
DoubleMatrix2D omega = sigmaInv.assign(
sigmaInv.zMult(s, null).zMult(sigmaInv, null),
minus
);
DoubleMatrix2D omegaLambda = omega.zMult(lambda, null);
// lambda梯度
DoubleMatrix2D dLambda = omegaLambda.zMult(phi, null);
// phi梯度
DoubleMatrix2D dPhi = lambda.viewDice().zMult(omegaLambda, null);
lambda.assign(dLambda, (x, y) -> x == 0 ? 0 : x - y * stepSize);
phi.assign(dPhi, (x, y) -> x == 1 ? 1 : x - y * stepSize);
// error的梯度直接为omega
error.assign(omega, (x, y) -> x == 0 ? 0 : x - y * stepSize);
}
广义混合效应模型（多层广义线性模型）之MCMC参数估计
这里的的广义混合效应模型是最简单的广义混合效应模型，关于混合效应模型，可以参考作者专栏写的另一篇文章
换一种写法
是固定截距，
是随机截距
我们通过MCMC吉布斯采样进行参数估计，整体过程如下
，链长，burn，thin，设
，
（2），
的设为0向量
（3）随机采样，
（4）计算后验似然函数与
后验似然函数的比值
（5）随机采样，
（6） 若，
（7）随机采样，
（8）计算后验似然函数与
后验似然函数的比值
（9）随机采样，
（10） 若，
（11）重复（3）到（10）直至迭代结束
部分代码如下
首先，我们需要实现一个logistic函数
/**
* logistic函数
* p(z) = e^z / (1 + e^z)
* @param intercept 固定截距
* @param z 随机截距
* @return logistic值矩阵
*/
public DoubleMatrix2D logistic(DoubleMatrix1D intercept, DoubleMatrix1D z) {
DoubleMatrix2D p = plusOuter(z, intercept);
p.assign(a -> {
double val = 1d / (1d + Math.exp(-a));
if (val <= 0)
return 1e-10;
else if (val >= 1d)
return 1 - 1e-10;
return val;
});
return p;
}
注意到一个小技巧，为了防止数值溢出，作者给logistic值设定了1个上界和下界，这是数值分析常用的做法。
colt 并没有实现不同形状向量相加方法，导致我们不得不写一个类似于外积的外和方法，主要是求
，实现如下
/**
* 类似于外积的外和 z[i, j] = x[i] + y[j]
*/
private DoubleMatrix2D plusOuter(DoubleMatrix1D x, DoubleMatrix1D y) {
int rows = x.size();
int columns = y.size();
DoubleMatrix2D out = x.like2D(rows,columns);
for (int row = rows; --row >= 0; ) out.viewRow(row).assign(y);
for (int column = columns; --column >= 0; ) out.viewColumn(column).assign(x, Functions.plus);
return out;
}
完整的实现和测试代码地址如下
inuyasha11/statsgitee.com
