BFGS算法的JAVA实现，，部分代码！（修正）

L-BFGS 主要思想，通过记忆m次的梯度信息来重构第k+1次的Hessian矩阵，用cpu的计算时间换取内存空间

package Regression;import java.io.BufferedWriter;import java.io.FileReader;import java.io.File;import java.io.IOException;import java.text.*;import java.util.*;import java.math.*;import weka.core.Instance;import weka.core.Instances;import weka.filters.Filter;import weka.filters.unsupervised.attribute.Remove;import weka.core.matrix.*;public class L_BFGS_Regression {static InstancesProcess inp=new InstancesProcess();Matrix theta=null;Matrix x=null;Matrix y=null;double batch_learningrate=0.00003;double lambda=0;double Max_Iter_Num=2000;double Cost_threashold_Rate=0.0001;double Cost_threashold=1;Matrix L_BFGS_Hessian=null;Matrix L_BFGS_y=null;Matrix L_BFGS_s=null;int L_BFGS_k=0;double[] lou_vector=null;// y=theta' * x/** * 返回损失函数 * @author 杭州趋数网络科技有限公司数据分析/数据挖掘工程师花京华 * */public double costfunction(){double costerr=0;int[] index_set={0};for(int i=0;i<x.getRowDimension();i++){double err=linearinstanceerror(i);costerr+=err*err; }costerr+=theta.normF();costerr=costerr/(2*x.getRowDimension());return costerr;}/** * 梯度向量 * 对每个样本执行Tao(J)/Tao(theta(j)) ，对线性回归 =(theta * x-y)* xj * @author 杭州趋数网络科技有限公司数据分析/数据挖掘工程师花京华 * * */public Matrix gradient(int i){Matrix grandient_delta=new Matrix(theta.getColumnDimension(),1,0);for(int j=0;j<theta.getColumnDimension();j++){int[] index_set={j};double Instanceserror=linearinstanceerror(i)*x.get(i, j)+lambda*theta.get(0, j);grandient_delta.set(j, 0, Instanceserror);}return grandient_delta;}public void linearLBFGS_descend_Batch(int m){L_BFGS_y=new Matrix(theta.getColumnDimension(),m,0);L_BFGS_s=new Matrix(theta.getColumnDimension(),m,0);lou_vector=new double[m];for(int i=0;i<m;i++)lou_vector[i]=0;L_BFGS_Hessian=Matrix.identity(theta.getColumnDimension(),theta.getColumnDimension());for(int it=0;it<Max_Iter_Num;it++){Matrix g1=gradient();//inp.display_matrix(g1);Matrix dk=L_BFGS_Hessian.times(g1).times(batch_learningrate);//dk is xk+1- xktheta.minusEquals(dk.transpose());//theat k+1Matrix g2=gradient();Matrix delta_gradient=g2.minus(g1);L_BFGS_k++;Hessian(dk,delta_gradient,m);//inp.display_matrix(delta_gradient);//inp.display_matrix(dk);//inp.display_matrix(L_BFGS_Hessian);inp.display_matrix(theta);double errcost=costfunction();System.out.println("ErrCost:"+errcost);}}/** * v'*M*v= V为N*m为矩阵 * */public Matrix vtMvMatrix(Matrix v,Matrix M){Matrix Middle=v.transpose().times(M).times(v);return Middle;}/** * L_BFGS 的头项 * @author 杭州趋数网络科技有限公司数据分析/数据挖掘工程师花京华 * @param y N*m维的梯度差分矩阵 * @param s N*m维的的待求解向量差分矩阵 * @param n 记忆次数 * */public Matrix Hessian_Header(Matrix y,Matrix s,int n){int[] index_set={0};int m=0;if(L_BFGS_k<=n-1){m=L_BFGS_k;}else{m=n;}Matrix I= Matrix.identity(y.getRowDimension(), y.getRowDimension());Matrix vector_y=y.getMatrix(0, y.getRowDimension()-1, index_set);Matrix vector_s=s.getMatrix(0, s.getRowDimension()-1, index_set);Matrix v=I.minus(vector_s.times(vector_y.transpose()).times(lou_vector[0]));Matrix Header=vtMvMatrix(v,I);//v0'Hv0for(int i=1;i<m;i++){index_set[0]=i;vector_y=y.getMatrix(0, y.getRowDimension()-1, index_set);vector_s=s.getMatrix(0, s.getRowDimension()-1, index_set);v=I.minus(vector_s.times(vector_y.transpose()).times(lou_vector[i]));Header=vtMvMatrix(v,Header);}return Header;}/** * L_BFGS 的尾项 * @author 杭州趋数网络科技有限公司数据分析/数据挖掘工程师花京华 * @param y N*m维的梯度差分矩阵 * @param s N*m维的的待求解向量差分矩阵 * @param n 记忆次数 * */public Matrix Hessian_Tail(Matrix y,Matrix s,int n){int m=0;if(L_BFGS_k<=n-1){m=L_BFGS_k;}else{m=n;}int[] index_set={m-1};Matrix I= Matrix.identity(y.getRowDimension(), y.getRowDimension());Matrix vector_s=s.getMatrix(0, s.getRowDimension()-1, index_set);Matrix Tail=vector_s.times(vector_s.transpose()).times(lou_vector[index_set[0]]);for(int j=0;j<m-2;j++)//for each item ,s0,…..sk-1 as middle {index_set[0]=j;vector_s=s.getMatrix(0, s.getRowDimension()-1, index_set);Matrix Middle=vector_s.times(vector_s.transpose()).times(lou_vector[j]);for(int i=j+1;i<m-2;i++)//for each v,v1….vk{index_set[0]=i;Matrix vector_v_s=s.getMatrix(0, s.getRowDimension()-1, index_set);Matrix vector_v_y=y.getMatrix(0, y.getRowDimension()-1, index_set);Matrix v=I.minus(vector_v_s.times(vector_v_y.transpose()).times(lou_vector[0]));Middle=vtMvMatrix(v,Middle);// ∏ v'*middle*v}Tail.plus(Middle);//∑ each item v(k)'*….*v(1)'*middle*v(1)….*v(k)}return Tail;}/** * L_BFGS 主要迭代开始 * @author 杭州趋数网络科技有限公司数据分析/数据挖掘工程师花京华 * @param vector_y 第k次的梯度差分 * @param vector_s 第k次的待求解向量差分 * @param n 记忆次数 * */public void Hessian(Matrix vector_y,Matrix vector_s,int n){refresh_s_y(vector_y,vector_s);refresh_lou_vector(vector_y,vector_s);L_BFGS_Hessian=Hessian_Header(L_BFGS_y,L_BFGS_s,n).plus(Hessian_Tail(L_BFGS_y,L_BFGS_s,n));}public static void main(String[] args) throws Exception {// TODO Auto-generated method stubString filename="E://拟合测试.csv";Instances org=inp.getInstances_fromcsv(filename);org.setClassIndex(org.numAttributes()-1);Gradient gd=new Gradient();gd.SetInput(org);//gd.lineargradientdescend_Batch();gd.linearLBFGS_descend_Batch(10);inp.display_matrix(gd.theta);}}

/********************************************************/

杭州趋数网络科技有限公司数据分析/数据挖掘工程师花京华

每一件事与每一个美丽景色，都有可能成为一生中的难忘。

相关文章：

你感兴趣的文章：

标签云：