矩阵常用函数的实现:
1、Java
1.1、函数的定义
package jm.app.algebra;
import org.ojalgo.matrix.PrimitiveMatrix;
import org.ojalgo.matrix.decomposition.Eigenvalue;
import java.math.BigDecimal;
import java.util.*;
public class AlgebraUtil {
//常用函数
// 1、协方差函数
public final static BigDecimal covariance(Matrix x, Matrix y) {
if (x == null || y == null || x.getRowNum() != y.getRowNum()) {
return null;
}
// 步骤1:计算X和Y的均值
int n = x.getRowNum();
BigDecimal xMean = new BigDecimal(0.0);
for (int i = 0; i < n; i++) {
xMean = xMean.add(x.getValue(i, 0));
}
xMean = xMean.multiply(new BigDecimal(1.0 / n));
BigDecimal yMean = new BigDecimal(0.0);
for (int i = 0; i < n; i++) {
yMean = yMean.add(y.getValue(i, 0));
}
yMean = yMean.multiply(new BigDecimal(1.0 / n));
// 步骤2:计算协方差矩阵
BigDecimal sum = new BigDecimal(0.0);
for (int i = 0; i < n; i++) {
sum = sum.add(x.getValue(i, 0).subtract(xMean).multiply(y.getValue(i, 0).subtract(yMean)));
}
BigDecimal value = sum.multiply(new BigDecimal(1.0 / (n - 1)));
return value;
}
//2、均值函数
public final static Matrix mean(Matrix x, int direction) {
if (x == null) {
return null;
}
if (direction == 0) {
// 新建矩阵newMatrix,用户保存计算后的结果
Matrix newMatrix = new Matrix(1, x.getColNum());
// 计算每一列的均值,并将其保存到newMatrix中
for (int c = 0; c < x.getColNum(); c++) {
BigDecimal mean = new BigDecimal(0.0);
for (int i = 0; i < x.getRowNum(); i++) {
mean = mean.add(x.getValue(i, c));
}
mean = mean.multiply(new BigDecimal(1.0 / x.getRowNum()));
newMatrix.setValue(0, c, mean);
}
return newMatrix;
} else if (direction == 1) {
// 新建矩阵newMatrix,用户保存计算后的结果
Matrix newMatrix = new Matrix(x.getRowNum(), 1);
// 计算每一列的均值,并将其保存到newMatrix中
for (int r = 0; r < x.getRowNum(); r++) {
BigDecimal mean = new BigDecimal(0.0);
for (int i = 0; i < x.getColNum(); i++) {
mean = mean.add(x.getValue(r, i));
}
mean = mean.multiply(new BigDecimal(1.0 / x.getColNum()));
newMatrix.setValue(r, 0, mean);
}
return newMatrix;
} else {
return null;
}
}
//3、最大值函数
public final static Matrix max(Matrix x, int direction) {
if (direction == 0) {
Matrix maxMat = new Matrix(1, x.getColNum());
for (int c = 0; c < x.getColNum(); c++) {
// 求出第c列最大值
BigDecimal maxValue = x.getValue(0, c);
for (int r = 0; r < x.getRowNum(); r++) {
BigDecimal value = x.getValue(r, c);
if (value.compareTo(maxValue) > 0) {
maxValue = value;
}
}
maxMat.setValue(0, c, maxValue);
}
return maxMat;
} else if (direction == 1) {
Matrix maxMat = new Matrix(x.getRowNum(), 1);
for (int r = 0; r < x.getRowNum(); r++) {
// 求出第c列最大值
BigDecimal maxValue = x.getValue(r, 0);
for (int c = 0; c < x.getColNum(); c++) {
BigDecimal value = x.getValue(r, c);
if (value.compareTo(maxValue) > 0) {
maxValue = value;
}
}
maxMat.setValue(r, 0, maxValue);
}
return maxMat;
} else {
return null;
}
}
//4、最小值函数
public final static Matrix min(Matrix x, int direction) {
if (direction == 0) {
Matrix minMat = new Matrix(1, x.getColNum());
for (int c = 0; c < x.getColNum(); c++) {
// 求出第c列最大值
BigDecimal minValue = x.getValue(0, c);
for (int r = 0; r < x.getRowNum(); r++) {
BigDecimal value = x.getValue(r, c);
if (value.compareTo(minValue) < 0) {
minValue = value;
}
}
minMat.setValue(0, c, minValue);
}
return minMat;
} else if (direction == 1) {
Matrix minMat = new Matrix(x.getRowNum(), 1);
for (int r = 0; r < x.getRowNum(); r++) {
// 求出第c列最大值
BigDecimal minValue = x.getValue(r, 0);
for (int c = 0; c < x.getColNum(); c++) {
BigDecimal value = x.getValue(r, c);
if (value.compareTo(minValue) < 0) {
minValue = value;
}
}
minMat.setValue(r, 0, minValue);
}
return minMat;
} else {
return null;
}
}
//5、归一化函数
public final static Matrix normalize(Matrix x, int direction) {
if (direction == 0) {
// 调用最大值函数和最小值函数得到maxmat和minmat
// maxmat表示x每一列的最大值
// maxmin表示x每一列的最小值
Matrix maxMat = AlgebraUtil.max(x, 0);
Matrix minMat = AlgebraUtil.min(x, 0);
// 开始对x的每一列进行归一化处理
Matrix normalizedMat = new Matrix(x.getRowNum(), x.getColNum());
for (int c = 0; c < x.getColNum(); c++) {
// 获取x的第c列的列向量
Matrix columnVector = getColumnVector(normalizeColumn(x, c, maxMat.getValue(0, c), minMat.getValue(0, c)), c);
normalizedMat = setColumnVector(normalizedMat, c, columnVector);
}
return normalizedMat;
} else if (direction == 1) {
Matrix maxMat = AlgebraUtil.max(x, 1);
Matrix minMat = AlgebraUtil.min(x, 1);
Matrix normalizedMat = new Matrix(x.getRowNum(), x.getColNum());
for (int r = 0; r < x.getRowNum(); r++) {
Matrix rowVector = getRowVector(normalizeRow(x, r, maxMat.getValue(r, 0), minMat.getValue(r, 0)), r);
normalizedMat = setRowVector(normalizedMat, r, rowVector);
}
return normalizedMat;
} else {
return null;
}
}
// 6、行列式函数
// 获取ojalog库中用于矩阵操作的工厂类
private final static PrimitiveMatrix.Factory matrixFactory = PrimitiveMatrix.FACTORY;
// 将ojalgo库中的矩阵类型转换为我们定义的Matrix类型
private final static Matrix convert(PrimitiveMatrix primitiveMatrix) {
Long rowNum = primitiveMatrix.countRows();
Long colNum = primitiveMatrix.countColumns();
Matrix matrix = new Matrix(rowNum.intValue(), colNum.intValue());
for (int i = 0; i < rowNum; i++) {
for (int j = 0; j < colNum; j++) {
matrix.setValue(i, j, primitiveMatrix.get(i, j));
}
}
return matrix;
}
// 将我们定义的Matrix类型转换为ojalgo库中的矩阵类型
private final static PrimitiveMatrix convert(Matrix matrix) {
int rowNum = matrix.getRowNum();
int colNum = matrix.getColNum();
PrimitiveMatrix primitiveMatrix = matrixFactory.makeZero(rowNum, colNum);
PrimitiveMatrix.DenseReceiver matrixBuilder = primitiveMatrix.copy();
for (int i = 0; i < rowNum; i++) {
for (int j = 0; j < colNum; j++) {
matrixBuilder.add(i, j, matrix.getValue(i, j).doubleValue());
}
}
return matrixBuilder.build();
}
public final static BigDecimal determinant(Matrix mat) {
if (mat.getColNum() != mat.getRowNum()) {
return null;
}
// 将mat矩阵转换为ojalgo库所使用的矩阵类型
PrimitiveMatrix primitiveMatrix = convert(mat);
// 直接使用ojalgo库的求行列式方法
return primitiveMatrix.getDeterminant().toBigDecimal();
}
//7、矩阵求逆函数
public final static Matrix inverse(Matrix mat) {
if (mat == null || mat.getRowNum() != mat.getColNum()) {
return null;
}
PrimitiveMatrix primitiveMatrix = convert(mat);
return convert(primitiveMatrix.invert());
}
//8、矩阵特征值和特征向量函数
public final static Matrix[] eigen(Matrix a) {
Matrix eigenvalues = new Matrix(a.getRowNum(), a.getColNum());
Matrix eigenvectors = new Matrix(a.getRowNum(), a.getColNum());
PrimitiveMatrix primitiveMatrix = convert(a);
List<Eigenvalue.Eigenpair> eigenpairList = primitiveMatrix.getEigenpairs();
for (int i = eigenpairList.size()-1; i >= 0; i--) {
// 按从小到大的顺序排列提取特征值
Eigenvalue.Eigenpair eigenpair = eigenpairList.get(i);
double eigenvalue = eigenpair.value.getReal();
int index = eigenpairList.size() - 1 - i;
eigenvalues.setValue(index, index, eigenvalue);
// 提取特征向量
for (int j = 0; j < eigenpair.vector.size(); j++) {
double element = eigenpair.vector.get(j).getReal();
eigenvectors.setValue(j, index, element);
}
}
Matrix[] valuesAndVectors = new Matrix[2];
valuesAndVectors[0] = eigenvalues;
valuesAndVectors[1] = eigenvectors;
return valuesAndVectors;
}
// 9、矩阵正交化函数
// 9.1 向量单位化
public final static Matrix unitize(Matrix x) {
if (x.getColNum() == 1) {
Matrix unitizedVector = new Matrix(x.getRowNum(), x.getColNum());
int elementNum = x.getRowNum();
BigDecimal denominator = l2Norm(x);
for (int i = 0; i < elementNum; i++) {
BigDecimal value = x.getValue(i, 0);
BigDecimal unitizedValue = value.multiply(new BigDecimal(1.0 / denominator.doubleValue()));
unitizedVector.setValue(i, 0, unitizedValue);
}
return unitizedVector;
} else if (x.getRowNum() == 1) {
Matrix unitizedVector = new Matrix(x.getRowNum(), x.getColNum());
int elementNum = x.getRowNum();
BigDecimal denominator = l2Norm(x);
for (int i = 0; i < elementNum; i++) {
BigDecimal value = x.getValue(0, i);
BigDecimal unitzedValue = value.multiply(new BigDecimal(1.0 / denominator.doubleValue()));
unitizedVector.setValue(i, 0, unitzedValue);
}
return unitizedVector;
} else {
return null;
}
}
// l2Norm方法实现响亮的第二范数,这里用于计算||x||
public final static BigDecimal l2Norm(Matrix x) {
if (x.getColNum() == 1) {
int elementNum = x.getRowNum();
BigDecimal sum = new BigDecimal(0.0);
for (int i = 0; i < elementNum; i++) {
BigDecimal value = x.getValue(i, 0);
sum = sum.add(value.multiply(value));
}
BigDecimal l2Value = new BigDecimal(Math.sqrt(sum.doubleValue()));
return l2Value;
} else if (x.getRowNum() == 1) {
int elementNum = x.getColNum();
BigDecimal sum = new BigDecimal(0.0);
for (int i = 0; i < elementNum; i++) {
BigDecimal value = x.getValue(0, i);
sum = sum.add(value.multiply(value));
}
BigDecimal l2Value = new BigDecimal(Math.sqrt(sum.doubleValue()));
return l2Value;
} else {
return null;
}
}
//9.2矩阵正交化
public final static Matrix orthogonalize(Matrix a) {
int dimension = a.getRowNum();
int vectorNum = a.getColNum();
Matrix b = new Matrix(dimension, vectorNum);
for (int i = 0; i < vectorNum; i++) {
Matrix ai = getColumnVector(a, i);
Matrix bi = copy(ai);
for (int j = 0; j < i; j++) {
Matrix bj = getColumnVector(b, j);
BigDecimal coefficient = inner(bj, ai).multiply(new BigDecimal(1.0 / inner(bj, bj).doubleValue()));
bi = substract(bi, dot(bj, coefficient));
}
// 将bi设置为矩阵b的第i列的列向量
b = setColumnVector(b, i, bi);
}
// 最后对矩阵b的每一列向量进行单位化处理
for (int i = 0; i < vectorNum; i++) {
b = setColumnVector(b, i, unitize(getColumnVector(b, i)));
}
return b;
}
}
1.2、函数的实现
package jm.app.algebra;
public class test {
public static void main(String[] args) {
// 1、设A=[6, 7, 1; 2, 2, 4],尝试利用本章给出的最大值函数和最小值函数分别计算A每行的最大值和最小值。
Matrix A = new Matrix(2, 3);
A.setValue(0,0,6);
A.setValue(0,1,7);
A.setValue(0,2,1);
A.setValue(1,0,2);
A.setValue(1,1,2);
A.setValue(1,2,4);
Matrix max = AlgebraUtil.max(A,1);
Matrix min = AlgebraUtil.min(A,1);
System.out.println(max);
System.out.println(min);
// 2、设A=[5, 2, 0; 2, 8, 3; 4 9 7],尝试利用本章给出的矩阵特征值和特征向量函数来计算A的所有特征值及其对应的特征向量。
Matrix B = new Matrix(3, 3);
B.setValue(0,0,5);
B.setValue(0,1,2);
B.setValue(0,2,0);
B.setValue(1,0,2);
B.setValue(1,1,8);
B.setValue(1,2,3);
B.setValue(2,0,4);
B.setValue(2,1,9);
B.setValue(2,2,7);
System.out.println(AlgebraUtil.eigen(B)[0]);
System.out.println(AlgebraUtil.eigen(B)[1]);
// 3、设A=[5, 2, 0; 2, 8, 3; 4 9 7],尝试利用本章给出的行列式函数求解出A的行列式。
System.out.println(AlgebraUtil.determinant(B));
// 4、设A=[5, 2, 0; 2, 8, 3; 4 9 7],尝试利用本章给出的矩阵求逆的实现来求解出A的逆矩阵A-1。
System.out.println(AlgebraUtil.inverse(B));
}
}
1.3、输出结果
7.000000 4.000000
1.000000 2.000000
13.260804 0.000000 0.000000 0.000000 4.218923 0.000000 0.000000 0.000000 2.520274
-0.130532 0.919425 0.606120 -0.539151 -0.359071 -0.751505 -0.858435 -0.160392 0.968602
141
0.205674 -0.099291 0.042553 -0.014184 0.248227 -0.106383 -0.099291 -0.262411 0.255319
Process finished with exit code 0
2、Python
2.1、函数的定义
import numpy as np
# 1、协方差函数
def covariance(a,b):
return np.cov(a,b)
# 2、均值函数
def mean(a,direction):
return np.mean(a,direction)
# 3、最大值函数
def max(a,direction):
return np.max(a,direction)
# 4、最小值函数
def min(a,direction):
return np.min(a,direction)
# 5、归一化函数
def normalize(a):
mx = max(a)
mn = min(a)
return [(float(i)-mn)/(mx-mn) for i in a]
# 6、行列式函数
def determinant(a):
return np.linalg.det(a)
# 7、矩阵求逆
def inverse(a):
return np.linalg.inv(a);
# 8、矩阵特征值和特征向量函数
def eigen(a):
eigenvalue,engenvector = np.linalg.eig(a)
# eigenvalue 按特征值由大到小排列
# engenvector 的第i列对应eigenvalue第i个特征值向量
return eigenvalue,engenvector
# 9、矩阵正交化函数
# 9.1 向量单位化
def unitize(a):
unitized_a = np.copy(a)
# 判断a是列向量还是行向量
if len(a.shape) == 1: #这种情况表示a为list类型
# 计算分母||a||
denominator = np.linalg.norm(a)
# 计算a/||a||
for i in range(len(a)):
unitized_a[i] = a[i] / denominator
elif a.shape[0] == 1: #这种情况表示a为行向量
# 计算分母||a||
denominator = np.linalg.norm(a)
# 计算a/||a||
for i in range(a.shape[1]):
unitized_a[0][i] = a[0][i] / denominator
elif a.shape[1] == 1: #这种情况表示a为行向量
# 计算分母||a||
denominator = np.linalg.norm(a)
# 计算a/||a||
for i in range(a.shape[0]):
unitized_a[i][0] = [i][0] / denominator
else:
print("仅可以单位化行向量或列向量")
return unitized_a
# 9.2矩阵正交化
def orthogonalize(a):
dimension = a.shape[0]
vectorNum = a.shape[1]
# 创建矩阵b
b = np.random.rand(dimension,vectorNum)
for i in range(vectorNum):
# 计算出ai和bi
ai = a[:,i]
bi = np.copy(ai)
for j in range(i):
bj = b[:,j]
coefficient = np.multiply(np.inner(bj,ai),1.0/np.inner(bi,bj))
coefficient =coefficient[0][0]
# 设置bi为矩阵b的第i列的列向量
b[:,i] = bi[:,0]
# 最后对矩阵的每一列向量进行单位化
for i in range(vectorNum):
b[:,i] = unitize(b[:,i])
return b
2.2、函数的实现
if __name__ == "__main__":
a = np.array([[6,7,1],[2,2,4]])
# 1、设A=[6, 7, 1; 2, 2, 4],尝试利用本章给出的最大值函数和最小值函数分别计算A每行的最大值和最小值。
print(max(a,1))
print(min(a,1))
# 2、设A=[5, 2, 0; 2, 8, 3; 4 9 7],尝试利用本章给出的矩阵特征值和特征向量函数来计算A的所有特征值及其对应的特征向量。
A = np.array([[5,2,0],[2,8,3],[4,9,7]])
print(eigen(A))
# 3.设A=[5, 2, 0; 2, 8, 3; 4 9 7],尝试利用本章给出的行列式函数求解出A的行列式。
print(determinant(A))
# 4.设A=[5, 2, 0; 2, 8, 3; 4 9 7],尝试利用本章给出的矩阵求逆的实现来求解出A的逆矩阵A-1。
print(inverse(A))
2.3、输出结果
[7 4]
[1 2]
(array([13.26080355, 4.21892288, 2.52027357]), array([[-0.1277134 , -0.91942511, 0.4431994 ],
[-0.52750766, 0.35907096, -0.54950663],
[-0.8398958 , 0.16039174, 0.70824907]]))
141.00000000000003
[[ 0.20567376 -0.09929078 0.04255319]
[-0.0141844 0.24822695 -0.10638298]
[-0.09929078 -0.26241135 0.25531915]]
Process finished with exit code 0
注:
此文中,用java语言求得的特征向量和Python求得的特征向量不一致,可能是Java语言求特征向量的程序写错了,但是目前没找到错误位置。