对了，这里面的奇异值分解方法我是直接抄袭与JAMA包。
package Matrix;
import java.io.Serializable;
import java.util.Arrays;
import MathUtil.MathUtil;
public final class Matrix implements Serializable {
/**
*
*/
private static final long serialVersionUID = -7914610542408796822L;
private int rownums;
private int colnums;
private Vector[] vectorArrays=null;
private Vector s=null;
private Matrix U=null;
private Matrix V=null;
// private double MatrixElements;
/***
* basic initialization...
* */
public Matrix(){}
/**
* initialization....
* */
public Matrix(int rownums,int colnums){
if(isInit(rownums,colnums)){
this.rownums=rownums;
this.colnums=colnums;
vectorArrays=new Vector[rownums];
init();
}else{
try{
throw new Exception("Error,please check this matrix");
}catch(Exception e){
e.printStackTrace();
}
}
}
/**
* initializing the value of specified position in this matrix...
* */
public double setMatrixValue(int rows,int cols,double MatrixElements){
if(isOutofBounds(rows,cols)){
try{
// throw new Exception("Error,please check the Matrix");
}catch(Exception e){
// e.printStackTrace();
}
}
return setValue(rows,cols,MatrixElements);
}
/**
* get the specified position of the matrix
* */
public double getMatrixValue(int rows,int cols){
if(isOutofBounds(rows,cols)){
try{
//throw new Exception("Error,please check the Matrix");
}catch(Exception e){
//e.printStackTrace();
}
}
return vectorArrays[rows-1].peek(cols-1);
}
/**
* construct the Matrix contains the value
* */
private double setValue(int rows,int cols,double MatrixElements){
return vectorArrays[rows-1].add(cols-1,MatrixElements);
}
/**
* the default value in each cell of this Matrix is 0;
* */
private Vector[] init(){
for(int i=1;i<=this.rownums;i++){
vectorArrays[i-1]=new Vector(colnums);
for(int j=1;j<=this.colnums;j++){
vectorArrays[i-1].add(j-1, 0);
}
}
return vectorArrays ;
}
/**
* toString
* */
public String toString(){
return Arrays.toString(vectorArrays);
}
private boolean isInit(int rows,int cols){
if(rows*cols<=0){
return false;
}else{
return true;
}
}
private boolean isOutofBounds(int rows,int cols){
if(rows<=0||rows>this.rownums&&cols<=0||cols>this.colnums){
return true;
}else{
return false;
}
}
/**
* transpose the matrix
* */
public Matrix tranposeMatrix(){
int rowsT=colnums;
int colsT=rownums;
Matrix matrixTranspose=new Matrix(rowsT,colsT);
for(int i=1;i<=rowsT;i++){
for(int j=1;j<=colsT;j++){
matrixTranspose.setMatrixValue(i, j, getMatrixValue(j,i));
}
}
return matrixTranspose;
}
/**
* the multiple of the matrix
* */
public Matrix mutiple(Matrix matrix){
int rows=matrix.getRownums();
int cols=matrix.getColnums();
Matrix matrixU=this.clone();
// System.out.println(matrixU.toString());
Matrix matrixA=null;
if(this.colnums==rows){
// matrixA.setRownums(rownums);
// matrixA.setColnums(cols);
matrixA=new Matrix(rownums,cols);
double d;
for(int i=1;i<=rownums;i++){
int j;
for(j=1;j<=cols;j++){
d=0;
for(int k=1;k<=colnums;k++){
d+=matrixU.getMatrixValue(i, k)*matrix.getMatrixValue(k, j);
// System.out.println("本矩阵的第"+i+"行第"+j+"列元素"+matrixU.getMatrixValue(i, k)+"乘被乘矩阵第"+j+"行第"+i+"列元素"+matrix.getMatrixValue(k, j)+"="+d);
// System.out.println("第i="+i+"次循环"+"第j="+j+"次循环");
// System.out.println("此时的d="+d);
// //if(j==2){
// //break;
// //}
}
matrixA.setMatrixValue(i, j, d);
}
}
}else{
try{
throw new Exception("计算失败，请确认该矩阵是否正确");
}catch(Exception e){
e.printStackTrace();
}
}
return matrixA;
}
// /**
// * two matrix multiply
// * */
// public Matrix mutiple(Matrix matrixA,Matrix matrixB){
// Matrix matrix=null;
// if(isMultiple(matrixA,matrixB)){
// matrix=new Matrix(matrixA.getRownums(),matrixB.getColnums());
// for(int i=1;i<=matrixA.getRownums();i++){
// double d=0;
// for(int j=1;j<=matrixA.getColnums();j++){
// d+=matrixA.getMatrixValue(i, j)*matrixB.getMatrixValue(j, i);
// }
// for(int j=1;j<=matrixB.getColnums();j++){
// matrix.setMatrixValue(i, j, d);
// }
// }
// }else{
// try{
// throw new Exception("计算失败，请确认矩阵是否正确");
// }catch(Exception e){
// e.printStackTrace();
// }
// }
// return matrix;
// }
/**
* the matrix equals
* @param precisions
* */
public boolean isEqual(Matrix matrix,double eps){
if(this.rownums==matrix.getRownums()&&this.colnums==matrix.getColnums()){
for(int i=1;i<=rownums;++i){
for(int j=1;j<=colnums;j++){
if((this.getMatrixValue(i, j)-matrix.getMatrixValue(i, j))<=eps){
return true;
}else{
return false;
}
}
}
}
return false;
}
/**
* the plus of matrix
*
* */
public Matrix plusMatrix(Matrix matrix){
Matrix matrixSum=new Matrix(this.rownums,this.colnums);
if(this.rownums==matrix.rownums&&this.colnums==matrix.colnums){
for(int i=1;i<=this.rownums;i++){
for(int j=1;j<=this.colnums;j++){
matrixSum.setMatrixValue(i, j, this.getMatrixValue(i, j)+matrix.getMatrixValue(i, j));
}
}
}else{
try{
throw new Exception("请注意，该矩阵无法相加");
}catch(Exception e){
e.printStackTrace();
}
}
return matrixSum;
}
public Matrix minusMatrix(Matrix matrix){
Matrix matrixminus=new Matrix(this.rownums,this.colnums);
if(this.rownums==matrix.rownums&&this.colnums==matrix.colnums){
for(int i=1;i<=this.rownums;i++){
for(int j=1;j<this.colnums;j++){
matrixminus.setMatrixValue(i, j, this.getMatrixValue(i, j)-matrix.getMatrixValue(i, j));
}
}
}else{
try{
throw new Exception("请注意，该矩阵无法相加");
}catch(Exception e){
e.printStackTrace();
}
}
return matrixminus;
}
// /***
// * @param eps the precisions
// * */
// public boolean isEqual(Matrix matrixA,Matrix matrixB,double eps){
// if(matrixA.isEqual(matrixB, eps)){
// return true;
// }else{
// return false;
// }
//
// }
/**
* return the boolean of two matrix,multiply by each matrix
* */
private boolean isMultiple(Matrix matrix){
int rows=matrix.getRownums();
//int cols=matrix.getColnums();
if(this.colnums==rows){
return true;
}else{
return false;
}
}
// private boolean isMultiple(Matrix matrixA,Matrix matrixB){
// if(matrixA.isMultiple(matrixB)){
// return true;
// }else{
// return false;
// }
// }
/**
* multiple the constant
* */
public Matrix multipleConstant(double K){
for(int i=1;i<=this.rownums;i++){
for(int j=1;j<=this.colnums;j++){
double d=this.getMatrixValue(i, j)*K;
this.setMatrixValue(i, j, d);
}
}
return this;
}
/**
* vector mutiple by the constant
* */
public Vector mutipleConstantV(double K,int cols){
if(cols>colnums||cols<=0){
try{
throw new Exception("请确认是否有误");
}catch(Exception e){
e.printStackTrace();
}
}
for(int i=1;i<=rownums;i++){
this.setMatrixValue(i, cols, K*this.getMatrixValue(i, cols));
}
return this.getVerticalVector(cols);
}
/**
* vector mutiple by the constant
* */
public Vector mutipleConstantH(double K,int rows){
if(rows>rownums||rows<=0){
try{
throw new Exception("请确认是否有误");
}catch(Exception e){
e.printStackTrace();
}
}
for(int i=1;i<=colnums;i++){
this.setMatrixValue(rows, i, K*this.getMatrixValue(rows, i));
}
return this.getHorizontalVector(rows);
}
/**
* get the egenValue of this Matrix
* how? Gaussian 消元法
* */
public double egenValue(){
double value=1;
if(this.rownums==this.colnums){
for(int j=1;j<=colnums;j++){
double max=this.getMaxAllRowValue(j);
if(Math.abs(max)!=0){
int rows=this.getMaxValueRow(max, j);
for(int i=1;i<this.getRownums();i++){
this.SwapMatrixValue(i, j, rows, j);
Vector vector=this.getHorizontalVector(i+1);
Vector vectorH=mutipleConstantH((1/max),i); //归一化
Vector vectorU=vector.vecMinus(vectorH.vecMutipleConstant(vectorH.peek(j-1)));
this.setMatrixHVector(i, vectorU);
}
}
}
for(int i=1;i<=rownums;i++){
value*=this.getMatrixValue(i, i);
}
}else{
try{
throw new Exception("不是行列式，请确认");
}catch(Exception e){
e.printStackTrace();
}
}
return value;
}
// /**
// * 高斯消元法
// * */
// public double egenDet(){
// double value=1;
// for(int j=1;j<=colnums;j++){
System.out.println("-----第"+j+"次行列式变化---"+this.toString());
// double Max_ABS=this.getABSMaxGivenRowValue(j);
System.out.println("-------这时的单列最大值------"+Max_ABS);
// int rows_Max_Abs=this.getMaxValueRow(Max_ABS, j);
// Vector vector=this.getHorizontalVector(rows_Max_Abs);
System.out.println("这时的最大值的行向量为"+vector.toString());
// if(Max_ABS==0){
// value=0;
// break;
// }else{
//
// for(int i=j;i<=rownums;i++){
// Vector unitVector=(vector.vecMutipleConstant(1/(Max_ABS)));
// System.out.println("此时的unitVector="+unitVector.toString());
/**
* 下面程序用来解决消元向量的绝对值问题
* */
// double unitVectorFirst=Math.abs(unitVector.peek(i-1));
// unitVector.add(i-1, unitVectorFirst);
// System.out.println("此时的unitVector="+unitVector.toString());
// System.out.println("--------------------------------------------追踪这时的行变化:"+i+"");
//
System.out.println("开始计算第"+j+"列unitVector="+unitVector.toString());
// Vector eachVector=this.getHorizontalVector(i);
// double Parameter=unitVector.peek(i-1);
System.out.println("此时的Parameter"+Parameter);
// double p=eachVector.peek(i-1);
System.out.println("此时的eachVector"+eachVector.toString()+"它的第"+i+"个向量值为"+eachVector.peek(i-1));
// if(i==rows_Max_Abs){
// continue;
// }else{
System.out.println("开始计算第"+j+"列"+"第"+i+"行"+eachVector.toString());
if(Parameter<0){
//
// Vector newUnitVector=unitVector.vecMutipleConstant(p);
System.out.println("乘上P以后新的消元向量为"+newUnitVector.toString());
// Vector neweachVector=null;
System.out.println("第"+j+"列最大值"+Max_ABS);
System.out.println("此时的newUnitVector="+newUnitVector.toString());
// neweachVector=eachVector.vecMinus(newUnitVector);
System.out.println("消元后的neweachVector="+neweachVector.toString());
System.out.println("归一化后的每个向量"+neweachVector.toString());
// this.setMatrixHVector(i, neweachVector);
System.out.println("-----第"+j+"次行列式变化---"+this.toString());
// this.SwapMatrixRows(j, rows_Max_Abs);
// this.multipleConstant(-1);
System.out.println("第"+i+"次变化时该行列式"+this.toString());
//
}else{
System.out.println("------------------------------该语句执行了啊-------------------------------");
Vector newUnitVector=unitVector.vecMutipleConstant(p);
System.out.println("乘上P以后新的消元向量为"+newUnitVector.toString());
Vector neweachVector=null;
System.out.println("第"+j+"列最大值"+Max_ABS);
System.out.println("此时的newUnitVector="+newUnitVector.toString());
neweachVector=eachVector.vecMinus(newUnitVector);
System.out.println("消元后的neweachVector="+neweachVector.toString());
// System.out.println("归一化后的每个向量"+neweachVector.toString());
this.setMatrixHVector(i, neweachVector);
System.out.println("-----第"+j+"次行列式变化---"+this.toString());
this.SwapMatrixRows(j, rows_Max_Abs);
this.multipleConstant(-1);
}
System.out.println("第"+i+"次变化时该行列式"+this.toString());
// }
// }
// }
System.out.println("第"+j+"次变化的行列式最终形式为"+this.toString());
// }
// for(int i=1;i<=rownums;i++){
// value*=this.getMatrixValue(i, i);
// }
System.out.println("最终的行列式"+this.toString());
// return value;
// }
// /**
// * 以下用来计算该矩阵的指定行最大值，并实现交换的过程方法
// * 假设这列最大值为aij
// * 假设实现到了rows行=cols列
// * rows=k
// * k<=i,j<=n;
// * 然后交换K行与i列
// * 交换K列与K行
// * */
// /**
// * 根据求出的最大值编写一个布尔类型
// * 来进行消元
// * 定义一个布尔型从s行到n行选择最大的元素作为主元的函数。
// * 并且交换
// * @param cols 不定行或者列
// * */
// private boolean findMaxofMatrix(int cols){
// if(this.getABSMaxGivenRowValue(cols)==0){
// return false;
// }else{
// this.getGaussianMatrixCols(cols);
// return true;
// }
// }
// /**
// * 消元计算生成上三角形
// * */
// private Matrix GaussianCalMatrix(int cols){
// int i;
// int j;
// int rows=cols;
// double parameter=this.getMatrixValue(rows, cols);
// Matrix unitMatrix=this.unitMatrix();
// for(i=rows;i<=rownums;i++){
// for(j=1;j<=colnums;j++){
// this.setMatrixValue(i, j, (this.getMatrixValue(i, j)-this.getMatrixValue(rows, j))*parameter);
// unitMatrix.setMatrixValue(i, j, (unitMatrix.getMatrixValue(i, j)-unitMatrix.getMatrixValue(rows,j))*parameter);
// }
// }
// return this;
// }
// /***
// * 以下方法构造一个单位行列式
// * */
public Matrix unitMatrix(){
Matrix matrix=new Matrix(rownums,colnums);
if(rownums==colnums){
for(int i=1;i<=rownums;i++){
System.out.println("i="+i+",");
for(int j=1;j<colnums;j++){
matrix.setMatrixValue(i, i, 1);
// System.out.println("i="+i+"j="+j);
}
}
}else{
try{
throw new Exception("不好意思，这不是个行列式的单位矩阵");
}catch(Exception e){
e.printStackTrace();
}
}
return matrix;
}
/**
* 矩阵的存储，典型的高斯消元法
* 以下程序设计思路，用两个数组用来存储所有的最大值的行或者列
* 三层嵌套for循环有点接近类似于之前的打印三角形，不太容易理解。
* */
public Matrix getAdverMatrix(){
int[] A=new int[this.rownums];
int[] B=new int[this.colnums];
int i,j,k,N;
double D;
if(this.rownums!=this.colnums){
try{
throw new Exception("请确认该矩阵是否为行列式");
}catch(Exception e){
e.printStackTrace();
}
}
Matrix matrix=this.clone();
N=rownums;
for(k=1;k<=N;k++){
D=0;
for(i=k;i<=N;i++){
for(j=k;j<=N;j++){
if(Math.abs(matrix.getMatrixValue(i, j))>D){
D=Math.abs(matrix.getMatrixValue(i, j));
A[k-1]=i;
B[k-1]=j;
}
}
}
if(D+1==1){
System.out.println("该矩阵不可逆");
break;
}
if(A[k-1]!=k){
for(j=1;j<=N;j++){
matrix.SwapMatrixValue(k, j, A[k-1], j);
}
}
if(B[k-1]!=k){
for(i=1;i<=N;i++){
matrix.SwapMatrixValue(i, k, i, B[k-1]);
}
}
double f=matrix.getMatrixValue(k, k);
matrix.setMatrixValue(k, k, (1/f));
for(j=1;j<=N;j++){
if(j!=k){
double d1=matrix.getMatrixValue(k, j);
double d2=matrix.getMatrixValue(k, k);
matrix.setMatrixValue(k, j, d1*d2);
}
}
for(i=1;i<=N;i++){
if(i!=k){
for(j=1;j<=N;j++){
if(j!=k){
double d=matrix.getMatrixValue(i, j);
double d1=matrix.getMatrixValue(i, k);
double d2=matrix.getMatrixValue(k, j);
matrix.setMatrixValue(i, j, d-d1*d2);
}
}
}
}
for(i=1;i<=N;i++){
if(i!=k){
double d1=matrix.getMatrixValue(i, k);
double d2=matrix.getMatrixValue(k, k);
matrix.setMatrixValue(i, k, -1*d1*d2);
}
}
}
for(k=N;k>=1;k--){
if(B[k-1]!=k){
for(j=1;j<=N;j++){
matrix.SwapMatrixValue(k, j, B[k-1], j);
}
}
if(A[k-1]!=k){
for(i=1;i<=N;i++){
matrix.SwapMatrixValue(i, k, i, A[k-1]);
}
}
}
return matrix;
}
/**
* 高斯消元法求行列式值
* */
public double egenDet(){
int i,j,k,nIs=1,nJs=1,N;
double f,det,q,d;
f=1;
det=1;
N=this.rownums;
if(N!=this.colnums){
try{
throw new Exception("请确认该矩阵是否是行列式");
}catch(Exception e){
e.printStackTrace();
}
}
/**
* 选取主元
* */
for(k=1;k<N;k++){
q=0;
for(i=k;i<=N;i++){
for(j=k;j<=N;j++){
d=Math.abs(this.getMatrixValue(i, j));
if(d>q){
q=d;
nIs=i;
nJs=j;
}
}
}
if(q+1==q){
det=0;
break;
}
if(nIs!=k){
f=-f;
for(j=k;j<=N;j++){
this.SwapMatrixValue(k, j, nIs, j);
}
}
if(nJs!=k){
f=-f;
for(i=k;i<=N;i++){
this.SwapMatrixValue(i, nJs, i, k);
}
}
det=det*this.getMatrixValue(k, k);
for(i=k+1;i<=N;i++){
double e=this.getMatrixValue(k, k);
d=this.getMatrixValue(i, k)/e;
for(j=k+1;j<=N;j++){
double d1=this.getMatrixValue(i, j);
double d2=this.getMatrixValue(k, j);
this.setMatrixValue(i, j, d1-d*d2);
}
}
det=f*det*this.getMatrixValue(rownums, colnums);
}
return det;
}
/**
* 以下求矩阵的特征值或者特征向量
* 用雅克比法计算矩阵的特征值和特征向量
* @param eps 计算精度
* @param nMaxItNum 允许迭代过程中最大迭代次数
* */
public boolean MJacobiEigenv(double eps,int nMaxItNum){
int i,j,p=1,q=1,l;
// nMaxItNum=1000;
double fm,cn=1,sn=1,omega,x,y,d;
l=1;
int state=-1;
int N=this.rownums;
if(this.rownums!=this.colnums){
try{
throw new Exception("不好意思，这个矩阵无法求特征值");
}catch(Exception e){
e.printStackTrace();
}
}
Matrix matrix=this.clone();
Matrix eigenVMatrix=new Matrix(N,N);
for(i=1;i<=N;i++){
eigenVMatrix.setMatrixValue(i, i, 1);
for(j=1;j<=N;j++){
if(i!=j){
eigenVMatrix.setMatrixValue(i, j, 0);
}
}
}
while(true){
fm=0;
for(i=2;i<=N;i++){
for(j=1;j<=i-1;j++){
d=Math.abs(matrix.getMatrixValue(i, j));
if(i!=j&&d>fm){
fm=d;
p=i;
q=j;
}
}
}
//这里可能需要些代码
System.out.println("state="+state);
if(fm<eps){
state=1;
return true;
// break;
}
if(l>nMaxItNum){
System.out.println("l="+l);
// System.out.println(eigenVMatrix.toString());
return false;
}
l=l+1;
x=-1*matrix.getMatrixValue(p, q);
y=(matrix.getMatrixValue(q, p)-matrix.getMatrixValue(p, q))/2;
omega=x/Math.sqrt(x*x+y*y);
if(y<0){
omega=-omega;
sn=1+Math.sqrt(1-omega*omega);
sn=omega/Math.sqrt(2*sn);
cn=Math.sqrt(1-sn*sn);
fm=matrix.getMatrixValue(p, p);
double d1=matrix.getMatrixValue(q, q);
double d2=matrix.getMatrixValue(p, q);
double d3=fm*cn*cn+d1*sn*sn+d2*omega;
matrix.setMatrixValue(q, q, d3);
double d4=matrix.getMatrixValue(q, q);
double d5=matrix.getMatrixValue(p, q);
double d6=fm*sn*sn+d4*cn*cn-d5*omega;
matrix.setMatrixValue(q, q, d6);
matrix.setMatrixValue(p, q, 0);
matrix.setMatrixValue(q, p, 0);
}
for(j=1;j<=N;j++){
if(j!=p&&j!=q){
fm=matrix.getMatrixValue(p, j);
double d1=matrix.getMatrixValue(q, j);
double d2=fm*cn+d1*sn;
matrix.setMatrixValue(p, j, d2);
double d3=matrix.getMatrixValue(q,j);
double d4=-1*sn+d3*cn;
matrix.setMatrixValue(q, j, d4);
}
}
for(i=1;i<=N;i++){
if(i!=p&&i!=q){
fm=matrix.getMatrixValue(i, p);
double d1=matrix.getMatrixValue(i, q);
double d2=fm*cn+d1*sn;
matrix.setMatrixValue(i, p, d2);
double d3=matrix.getMatrixValue(i, q);
double d4=-1*fm*sn+d3*cn;
matrix.setMatrixValue(i, q, d4);
}
}
for(i=1;i<=N;i++){
fm=eigenVMatrix.getMatrixValue(i, p);
double d1=eigenVMatrix.getMatrixValue(i, q);
double d2=fm*cn+d1*sn;
eigenVMatrix.setMatrixValue(i, p, d2);
double d3=-1*fm*sn+d1*cn;
eigenVMatrix.setMatrixValue(i, q, d3);
}
}
}
/**
* @param eps计算精度
* @param nMaxItNum 指定迭代次数
* */
public void getEigenV(double eps,int nMaxItNum){
boolean d;
d=this.MJacobiEigenv(eps,nMaxItNum);
if(d==false){
System.out.println(this.toString());
System.out.println("不好意思无法求解");
}else{
for(int i=1;i<=this.rownums;i++){
System.out.println(this.getMatrixValue(i, i));
System.out.println();
}
}
}
/**
* clone 方法
* */
public Matrix clone(){
Matrix matrix=new Matrix(this.rownums,this.colnums);
for(int i=1;i<=matrix.rownums;i++){
for(int j=1;j<=matrix.colnums;j++){
matrix.setMatrixValue(i, j, this.getMatrixValue(i, j));
}
}
return matrix;
}
/**
*
* 下面程序从cols列开始到列尾，求出最大绝对值所在行，和所在列
*
* */
public Matrix getGaussianMatrixCols(int cols){
int rows=cols;
double d=0;
d=this.getABSMaxGivenRowValue(cols);
for(int i=cols;i<=colnums;i++){
if(this.getABSMaxGivenRowValue(i)>d){
d=this.getABSMaxGivenRowValue(i);
}
}
int rowsofd=this.getMaxValueRow(d, cols);
int colsofd=0;
for(int j=cols;j<=colnums;j++){
if(Math.abs(this.getMatrixValue(rowsofd,j))==d){
colsofd=j;
}
}
Matrix unitMatrix=this.unitMatrix();
this.SwapMatrixRows(rows, rowsofd);
// this.SwapMatrixCols(cols, colsofd);
unitMatrix.SwapMatrixRows(rows,rowsofd);
return this;
}
/**
* 正交化矩阵方法
* */
public Matrix formularMatrix(){
Matrix matrix=this.clone();
Matrix matrixFormular=new Matrix(matrix.rownums,matrix.colnums);
if(matrix.egenDet()==0||rownums!=colnums){
try{
// throw new Exception("不好意思，无法正交化");
}catch(Exception e){
// e.printStackTrace();
}
}
// Matrix formularMatrix=new Matrix(matrix.rownums,matrix.colnums);
Vector Bvector1=matrix.getVerticalVector(1);
// formularMatrix.setMatrixVVector(1, Bvector1);
Vector A1=matrix.getVerticalVector(1);
Vector B1=A1;
double l=getParameter(B1,A1);
Vector B2=B1.vecMinus(B1.vecMutipleConstant(l));
matrixFormular.setMatrixVVector(1, B1.formVector());
matrixFormular.setMatrixVVector(2, B2.formVector());
for(int i=3;i<=matrix.getColnums();i++){
matrixFormular.setMatrixVVector(i, matrix.getSpecifiedFormularVector(i));
}
return matrixFormular;
}
/**
*
* @param Brpre上一个规范正交化向量
* @param Ar 待正交化的第r个列向量
* */
private double getParameter(Vector Brpre,Vector Ar){
double d=Brpre.vecMutipleVector(Ar);
double e=Brpre.vecMutipleVector(Brpre);
if(e==0){
try{
throw new Exception("除数出现0，请检查");
}catch(Exception e1){
e1.printStackTrace();
}
}
double f=-1*(d/e);
return f;
}
/**
*正交基过程
*@param Brpre
* @param Brcur当前规范正交化向量
* @param Brpre上一个规范正交化向量
* @param Ar 待正交化的第r个列向量
* */
private Vector getFormularVector(Vector Brpre,Vector Ar){
Vector vector=null;
double d=getParameter(Brpre,Ar);
vector=Brpre.vecMutipleConstant(d);
return vector;
}
/**
* 指定一个正交基向量的产生方法
* 采用递归方法或者for循环
* */
private Vector getSpecifiedFormularVector(int cols){
Matrix matrix=this.clone();
Vector Ar=this.getVerticalVector(cols);
Vector A1=this.getVerticalVector(1);
Vector B1=A1;
double l=getParameter(B1,A1);
Vector B2=B1.vecMinus(B1.vecMutipleConstant(l));
Vector temp;
for(int i=3;i<=cols;i++){
temp=B2;
double ll=getParameter(B2,matrix.getVerticalVector(i));
B2=B1.vecMinus(B2.vecMutipleConstant(ll));
B1=temp;
}
return B2;
}
/**
* return the max of the specified cols given the rows
*
* @param cols指定列的所有行的最大值
* */
public double getABSMaxGivenRowValue(int cols){
int rows=cols;
double max=Math.abs(this.getMatrixValue(rows, cols));
for(int i=rows+1;i<=this.rownums;i++){
if(Math.abs(this.getMatrixValue(i, cols))>max){
max=Math.abs(this.getMatrixValue(i, cols));
}
}
return max;
}
/**
* 获取指定列最大值的行值
* */
public int getMaxValueRow(double max,int cols){
int positions=1;
int rows=cols;
for(int i=rows;i<=this.rownums;i++){
if(Math.abs(this.getMatrixValue(i, cols))==max){
positions=i;
break;
}
}
return positions;
}
/**
*the following method is to Swap the values of the specified positions
* */
public Matrix SwapMatrixValue(int rows1,int cols1,int row2,int cols2){
double a=this.getMatrixValue(rows1, cols1);
double b=this.getMatrixValue(row2, cols2);
this.setMatrixValue(rows1, cols1, b);
this.setMatrixValue(row2, cols2, a);
return this;
}
/**
* to Swap the specified given rows between k,n
* @param k
* @param n
* */
public Matrix SwapMatrixRows(int k,int n){
Vector matrixk=this.getHorizontalVector(k);
Vector matrixn=this.getHorizontalVector(n);
this.setMatrixHVector(k, matrixn);
this.setMatrixHVector(n, matrixk);
return this;
}
/***
* to Swap the specified given cols between k,n
* */
public Matrix SwapMatrixCols(int k,int n){
Vector matrixk=this.getVerticalVector(k);
Vector matrixn=this.getVerticalVector(n);
this.setMatrixVVector(k, matrixn);
this.setMatrixVVector(n, matrixk);
return this;
}
/**
* 返回单一矩阵
* */
public Matrix getHorizontalMatrix(int rows){
Vector vector=getHorizontalVector(rows);
Matrix matrix=new Matrix(1,vector.length());
matrix.setMatrixHVector(1, vector);
return matrix;
}
/**
* the following method is to get the horizontal vectors;
* */
public Vector getHorizontalVector(int rows){
Vector vector=new Vector(this.colnums);
for(int i=1;i<=vector.length();i++){
vector.add(this.getMatrixValue(rows, i));
}
return vector;
}
public Vector getVerticalVector(int cols){
Vector vector=new Vector(this.rownums);
for(int i=1;i<=vector.length();i++){
vector.add(this.getMatrixValue(i, cols));
}
return vector;
}
/**
* 以下设置矩阵列向量
* */
public Matrix setMatrixVVector(int cols,Vector vector){
if(cols>colnums||vector.length()>rownums){
try{
throw new Exception("请确认该向量");
}catch(Exception e){
e.printStackTrace();
}
}
for(int i=0;i<vector.length();i++){
this.setMatrixValue(i+1, cols, vector.peek(i));
}
return this;
}
/**
* 以下设置矩阵行向量
* @param rows 指定行
* */
public Matrix setMatrixHVector(int rows,Vector vector){
if(rows>rownums||vector.length()>colnums){
try{
throw new Exception("请确认该向量");
}catch(Exception e){
e.printStackTrace();
}
}
for(int i=0;i<vector.length();i++){
this.setMatrixValue(rows, i+1, vector.peek(i));
}
return this;
}
/**
* 以下方法来自于JAMA包
* Singular Vector Decomposition
* @return
*/
public void getSVD(){
int m=this.rownums;
int n=this.colnums;
int nu = Math.min(this.rownums,this.colnums);
s=new Vector(Math.min(m+1, n));
U=new Matrix(m,nu);
V=new Matrix(n,n);
Matrix A=this.clone();
//Vector s = new Vector(Math.min(m+1,n));
int nct=Math.min(m-1, n);
int nrt=Math.max(0, Math.min(n-2,m));
Vector e=new Vector(n);
Vector work=new Vector(m);//构建一个长度为m的向量work
boolean wantu = true; //判断生成U
boolean wantv = true;
for(int k=0;k<Math.max(nct,nrt);k++){
if(k<nct){
s.add(k, 0);
for(int i=k;i<m;i++){
double d1=s.peek(k);
double d2=A.getMatrixValue(i+1, k+1);
double d3=MathUtil.getSquare(d1, d2);
s.add(k, d3);
}
if(s.peek(k)!=0){
if(A.getMatrixValue(k+1, k+1)<0){
double d1=s.peek(k);
double d2=-1*d1;
s.add(k,d2);
}
for (int i = k; i < m; i++) {
double d1=s.peek(k);
double d2=A.getMatrixValue(i+1, k+1);
A.setMatrixValue(i+1, k+1, d2/d1);
}
double dtmp=A.getMatrixValue(k+1, k+1);
A.setMatrixValue(k+1, k+1, dtmp+1);
}
double d1=s.peek(k);
s.add(k,-1*d1);
}
for (int j = k+1; j <n; j++) {
if ((k <nct) & (s.peek(k) != 0.0)) {
// Apply the transformation.
double t = 0;
for (int i = k; i <m; i++) {
t +=A.getMatrixValue(i+1, k+1)*A.getMatrixValue(i+1, j+1);
}
double d1=A.getMatrixValue(k+1, k+1);
t = -t/d1;
for (int i = k; i <m; i++) {
double dd=A.getMatrixValue(i+1, k+1);
double md=t*dd;
double dij=A.getMatrixValue(i+1, j+1);
double temp=dij+md;
A.setMatrixValue(i+1, j+1, temp);
}
}
// Place the k-th row of A into e for the
// subsequent calculation of the row transformation.
double tmp=A.getMatrixValue(k+1, j+1);
e.add(j, tmp);
}
if (wantu & (k <nct)) {
// Place the transformation in U for subsequent back
// multiplication.
for (int i = k; i < A.rownums; i++) {
double temp=A.getMatrixValue(i+1, k+1);
U.setMatrixValue(i+1, k+1, temp);
}
}
if (k < nrt) {
// Compute the k-th row transformation and place the
// k-th super-diagonal in e[k].
// Compute 2-norm without under/overflow.
e.add(k, 0);
for (int i = k+1; i < n; i++) {
double tmp1=e.peek(k);
double tmp2=e.peek(i);
double tmp=MathUtil.getSquare(tmp1, tmp2);
e.add(k,tmp);
}
if (e.peek(k) != 0.0) {
if (e.peek(k+1) < 0.0) {
// e[k] = -e[k];
double tmp1=e.peek(k);
double tmp2=-1*tmp1;
e.add(k, tmp2);
}
for (int i = k+1; i <n; i++) {
double tmp1=e.peek(k);
double tmpi=e.peek(i);
double tmp=tmpi/tmp1;
e.add(i, tmp);
}
double tmp=e.peek(k+1);
double tmp1=tmp+1;
e.add(k+1, tmp1);
}
double tmp=e.peek(k);
double tmpl=-1*tmp;
e.add(k, tmpl);
if ((k+1 <m) & (e.peek(k) != 0.0)) {
// Apply the transformation.
for (int i = k+1; i < m; i++) {
work.add(i, 0);
}
for (int j = k+1; j <n; j++) {
for (int i = k+1; i < m; i++) {
double tmp1=e.peek(j);
double tmp2=A.getMatrixValue(i+1, j+1);
double tmp3=work.peek(i);
double tmp4=tmp3+tmp1*tmp2;
work.add(i, tmp4);
}
}
for (int j = k+1; j < n; j++) {
double tmpej=e.peek(j)*(-1);
double tmpek1=e.peek(k+1);
double tmpn=tmpej/tmpek1;
double t = tmpn;
for (int i = k+1; i < m; i++) {
double tmpm=t*work.peek(i);
double tmpmm=A.getMatrixValue(i+1, j+1);
double tmppp=tmpmm+tmpm;
A.setMatrixValue(i+1, j+1, tmppp);
}
}
}
if (wantv) {
// Place the transformation in V for subsequent
// back multiplication.
for (int i = k+1; i < n; i++) {
double dd=e.peek(i);
V.setMatrixValue(i+1, k+1, dd);
}
}
}
}
// Set up the final bidiagonal matrix or order p.
int p = Math.min(n,m+1);
if (nct <n) {
double dtmp=A.getMatrixValue(nct+1, nct+1);
s.add(nct, dtmp);
}
if (m < p) {
s.add(p-1, 0);
}
if (nrt+1 < p) {
double dtmp=A.getMatrixValue(nrt+1, p);
e.add(nrt,dtmp);
}
e.add(p-1, 0);
// If required, generate U.
if (wantu) {
for (int j = nct; j < nu; j++) {
for (int i = 0; i < m; i++) {
U.setMatrixValue(i+1, j+1, 0);
}
U.setMatrixValue(j+1, j+1, 1);
}
for (int k = nct-1; k >= 0; k--) {
if (s.peek(k) != 0.0) {
for (int j = k+1; j < nu; j++) {
double t = 0;
for (int i = k; i < m; i++) {
double tmp1=U.getMatrixValue(i+1, k+1);
double tmp2=U.getMatrixValue(i+1, j+1);
double tmp3=tmp1*tmp2;
double tmp4=t+tmp3;
t=tmp4;
}
double UKK=U.getMatrixValue(k+1, k+1);
double mmm=-t/UKK;
t=mmm;
for (int i = k; i < m; i++) {
double Utmp=U.getMatrixValue(i+1, k+1);
double ttmp=t*Utmp;
double Uij=U.getMatrixValue(i+1, j+1);
double Utemp=Uij+ttmp;
U.setMatrixValue(i+1, j+1, Utemp);
}
}
for (int i = k; i < m; i++ ) {
double Utmp=U.getMatrixValue(i+1, k+1)*(-1);
U.setMatrixValue(i+1, k+1, Utmp);
}
double Ukktmp=U.getMatrixValue(k+1, k+1);
U.setMatrixValue(k+1, k+1, Ukktmp+1);
for (int i = 0; i < k-1; i++) {
U.setMatrixValue(i+1, k+1, 0);
}
} else {
for (int i = 0; i < m; i++) {
U.setMatrixValue(i+1, k+1, 0);
}
U.setMatrixValue(k+1, k+1, 1);
}
}
}
// If required, generate V.
//以下代码可能会出现问题
if (wantv) {
for (int k = n-1; k >= 0; k--) {
if ((k < nrt) & (e.peek(k) != 0.0)) {
for (int j = k+1; j < nu; j++) {
double t = 0;
for (int i = k+1; i <n; i++) {
double d1=V.getMatrixValue(i+1, k+1);
double d2=V.getMatrixValue(1+i, j+1);
double d3=d1*d2;
double d4=t+d3;
t=d4;
}
t = -t/V.getMatrixValue(k+2, k+1);
for (int i = k+1; i < n; i++) {
double d1=V.getMatrixValue(i+1, k+1);
double d2=V.getMatrixValue(i+1, j+1);
double d3=d1*t;
double d4=d3+d2;
V.setMatrixValue(i+1, j+1, d4);
}
}
}
for (int i = 0; i < n; i++) {
V.setMatrixValue(i+1, k+1, 0);
}
V.setMatrixValue(k+1, k+1, 1);
}
}
// Main iteration loop for the singular values. //下列程序迭代求解奇异值
int pp = p-1;
int iter = 0;
double eps = Math.pow(2.0,-52.0);
double tiny = Math.pow(2.0,-966.0);
while (p > 0) {
int k,kase;
// Here is where a test for too many iterations would go.
// This section of the program inspects for
// negligible elements in the s and e arrays. On
// completion the variables kase and k are set as follows.
// kase = 1 if s(p) and e[k-1] are negligible and k<p
// kase = 2 if s(k) is negligible and k<p
// kase = 3 if e[k-1] is negligible, k<p, and
// s(k), ..., s(p) are not negligible (qr step).
// kase = 4 if e(p-1) is negligible (convergence).
for (k = p-2; k >= -1; k--) {
if (k == -1) {
break;
}
if (Math.abs(e.peek(k)) <=
tiny + eps*(Math.abs(s.peek(k)) + Math.abs(s.peek(k+1)))) {
e.add(k, 0);
break;
}
}
if (k == p-2) {
kase = 4;
} else {
int ks;
for (ks = p-1; ks >= k; ks--) {
if (ks == k) {
break;
}
double t = (ks != p ? Math.abs(e.peek(ks)) : 0.) + //三目计算
(ks != k+1 ? Math.abs(e.peek(ks-1)) : 0.);
if (Math.abs(s.peek(ks)) <= tiny + eps*t) {
s.add(ks, 0);
break;
}
}
if (ks == k) {
kase = 3;
} else if (ks == p-1) {
kase = 1;
} else {
kase = 2;
k = ks;
}
}
k++;
// Perform the task indicated by kase.
switch (kase) {
// Deflate negligible s(p).
case 1: {
double f = e.peek(p-2);
e.add(p-2,0.0);
for (int j = p-2; j >= k; j--) {
double t = MathUtil.getSquare(s.peek(j),f);
double cs = s.peek(j)/t;
double sn = f/t;
s.add(j,t);
if (j != k) {
f = -sn*e.peek(j-1);
double etmp=e.peek(j-1);
double etmp2 = cs*etmp;
e.add(j-1, etmp2);
}
if (wantv) {
for (int i = 0; i < n; i++) {
t = cs*V.getMatrixValue(i+1, j+1) + sn*V.getMatrixValue(i+1, p);
double d1=V.getMatrixValue(i+1, j+1);
double d2=V.getMatrixValue(i+1, p);
double dtmp=-sn*d1+cs*d2;
V.setMatrixValue(i+1, p, dtmp);
V.setMatrixValue(i+1, j+1, t);
}
}
}
}
break;
// Split at negligible s(k).
case 2: {
double f = e.peek(k-1);
e.add(k-1,0.0);
for (int j = k; j < p; j++) {
double t = MathUtil.getSquare(s.peek(j),f);
double cs = s.peek(j)/t;
double sn = f/t;
s.add(j, t);
f = -sn*e.peek(j);
e.add(j ,(cs*e.peek(j)));
if (wantu) {
for (int i = 0; i < m; i++) {
double d1=U.getMatrixValue(i+1, j+1);
double d2=U.getMatrixValue(i+1, k);
t = cs*d1 + sn*d2;
U.setMatrixValue(i+1,k,-sn*d1 + cs*d2);
U.setMatrixValue(i+1, j+1, t);
}
}
}
}
break;
// Perform one qr step.
case 3: {
// Calculate the shift.
double scale = Math.max(Math.max(Math.max(Math.max(
Math.abs(s.peek(p-1)),Math.abs(s.peek(p-2))),Math.abs(e.peek(p-2))),
Math.abs(s.peek(k))),Math.abs(e.peek(k)));
double sp = s.peek(p-1)/scale;
double spm1 = s.peek(p-2)/scale;
double epm1 = e.peek(p-2)/scale;
double sk = s.peek(k)/scale;
double ek = e.peek(k)/scale;
double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
double c = (sp*epm1)*(sp*epm1);
double shift = 0.0;
if ((b != 0.0) | (c != 0.0)) {
shift = Math.sqrt(b*b + c);
if (b < 0.0) {
shift = -shift;
}
shift = c/(b + shift);
}
double f = (sk + sp)*(sk - sp) + shift;
double g = sk*ek;
// Chase zeros.
for (int j = k; j < p-1; j++) {
double t = MathUtil.getSquare(f,g);
double cs = f/t;
double sn = g/t;
if (j != k) {
e.add(j-1,t);
}
f = cs*s.peek(j) + sn*e.peek(j);
double dtmp=cs*e.peek(j) - sn*s.peek(j);
e.add(j, dtmp);
g = sn*s.peek(j+1);
double dtmp2=cs*s.peek(j+1);
s.add(j+1, dtmp2);
if (wantv) {
for (int i = 0; i < n; i++) {
t = cs*V.getMatrixValue(i+1, j+1) + sn*V.getMatrixValue(i+1, j+2);
double dtmpV=V.getMatrixValue(i+1, j+1)*(-sn);
double dtmpVi=V.getMatrixValue(i+1, j+2)*cs;
double dt=dtmpV+dtmpVi;
V.setMatrixValue(i+1, j+2, dt);
// V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
// V[i][j] = t;
V.setMatrixValue(i+1, j+1, t);
}
}
t = MathUtil.getSquare(f,g);
cs = f/t;
sn = g/t;
s.add(j,t);
f = cs*e.peek(j)+ sn*s.peek(j+1);
double dtmp1=cs*s.peek(j+1);
double dtmp3=(-sn)*e.peek(j);
double dtmp3t=dtmp1+dtmp3;
s.add(j+1,dtmp3t);
g = sn*e.peek(j+1);
double dddd = cs*e.peek(j+1);
e.add(j+1, dddd);
if (wantu && (j < m-1)) {
for (int i = 0; i < m; i++) {
double dtt=U.getMatrixValue(i+1, j+1);
double dttt=U.getMatrixValue(i+1, j+2);
t = cs*dtt + sn*dttt;
U.setMatrixValue(i+1,j+2,-sn*dtt + cs*dttt);
//U[i][j] = t;
U.setMatrixValue(i+1, j+1, t);
}
}
}
e.add(p-2, f);
iter = iter + 1;
}
break;
// Convergence.
case 4: {
// Make the singular values positive.
if (s.peek(k) <= 0.0) {
//数组转下标要小心
double d=(s.peek(k) < 0.0 ? -(s.peek(k)) : 0.0);
s.add(k, d);
if (wantv) {
for (int i = 0; i <= pp; i++) {
double dd=(-1)*V.getMatrixValue(i+1, k+1);
V.setMatrixValue(i+1, k+1, dd);
}
}
}
// Order the singular values.
while (k < pp) {
if (s.peek(k) >= s.peek(k+1)) {
break;
}
double t = s.peek(k);
double ddd=s.peek(k+1);
s.add(k, ddd);
s.add(k+1,t);
if (wantv && (k < n-1)) {
for (int i = 0; i < n; i++) {
V.SwapMatrixValue(i+1, k+2, i+1, k+1);
}
}
if (wantu && (k < m-1)) {
for (int i = 0; i < m; i++) {
V.SwapMatrixValue(i+1, k+2, i+1, k+1);
}
}
k++;
}
iter = 0;
p--;
}
break;
}
}
//return U;
}
/**
* 必须接在getSVD方法后使用
* */
public Matrix getU(){
return this.U;
}
/**
* 必须在getSVD方法后使用
* */
public Matrix getV(){
return this.V;
}
public Vector getSingularValues () {
return this.s;
}
/** Return the diagonal matrix of singular values
@return S
*/
public Matrix getS () {
Matrix X = new Matrix(this.colnums,this.colnums);
for (int i = 1; i <= X.rownums; i++) {
for (int j = 1; j <= X.colnums; j++) {
X.setMatrixValue(i, j, 0);
}
double tmp=s.peek(i-1);
X.setMatrixValue(i, i, tmp);
}
return X;
}
public int getRownums() {
return rownums;
}
public void setRownums(int rownums) {
this.rownums = rownums;
}
public int getColnums() {
return colnums;
}
public void setColnums(int colnums) {
this.colnums = colnums;
}
public static void main(String[] args){
Matrix matrix=new Matrix(2,2);
matrix.setMatrixValue(1, 1, 2);
matrix.setMatrixValue(2, 2, 2);
Matrix matrixA=new Matrix(2,2);
matrixA.setMatrixValue(1, 1, 20);
matrixA.setMatrixValue(2, 2, 20);
System.out.print(matrix.plusMatrix(matrixA).toString());
}
}这是我学习JavaSE的时候写的java版本矩阵类，提供了各种方法。
大学的线性代数很有用啊！！
1338
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
