参考参考了《Practical Optimization Methods With Mathematic Applications》中的8.4节中介绍的有效集法（Active-Set），有效集法只能优化中等规模的最优化问题，如果是大规模的最优化问题，应该采用其它算法，代码如下：（代码执行环境VS2010）
QP.h
#ifndef _QP_
#define _QP_
#pragma comment (lib,"./bin/blasd.lib")
#pragma comment (lib,"./bin/clapackd.lib")
#pragma comment (lib,"./bin/libf2cd.lib")
#pragma comment (lib,"./bin/tmglibd.lib")
#include <fstream>
#include <assert.h>
#include <vector>
#include <algorithm>
#include <iostream>
#include "Lapack.h"
using namespace std;
/*
Copy right
用拉格朗日乘子发求解等式约束的二次规划问题
求解问题为：min: 0.5*trans(x)*H*x + trans(c)*x; st: A*x = b
参考了《Practical Optimization Methods With Mathematic Applications》中的8.4节中介绍的有效集法（Active-Set）
*/
class QuadEquLag //等式约束问题
{
private:
//
double** rTmp;
public:
double* optimizer(double** H,int& m, int& n, double* c, double** A, double* b);//返回最优解
void calculateQ(double** reverseH, int& row, int& col, double** A,double** Q);
void calculateR(double** reverseH, int& row, int& col, double** A,double** R);
void twoMatrixMultiply(double** a, int& aRow, int& aCol, char* aType, double** b, int& bRow, int& bCol, char* bType,double** out);
void inverseMatrix(double** a, int& aRow, int& aCol); //计算矩阵a的逆矩阵
void testMatrixInverse();
};
class QuadProg : public QuadEquLag //不等式约束问题
{
private:
public:
//double* optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0);
double optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0); //有效集算法只能处理中等规模的数据
int rank(double** a, int& row, int& col, double tol=1e-10); //用SVD分解计算矩阵的秩
void solveLinearEquations(double** a, int& row, int& col, double*);
bool greaterThanZero(double* in, int& n);
bool equalZero(double* in, int& n);
void updateActiveSet(vector<int>& activeSet, double* mu);
double findAfak(vector<int>& activeSet, double* dk, double** A , double* b,double* xk, int& minIndex, int& m, int& n);
void updataXk(double* xk,double& afak, double* dk, int& n);
};
#endif
#include "QP.h"
void QuadEquLag::testMatrixInverse()
{
int N = 500;
double* A = new double[9];
int LDA = 500;
int* IPIV = new int[N];
int LWORK = N;
double* WORK = new double[N];
int INFO ;
ifstream readin;
readin.open("./data/m.txt",ios_base::in);
int k = 0;
if(readin.is_open())
{
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
double tmp = 0;
readin>>tmp;
A[k++] = tmp;
}
}
}
int m = N; int n = m;
dgetrf_(&m, &n, A, &m, IPIV, &INFO);
dgetri_(&N, A, &LDA, IPIV, WORK, &LWORK, &INFO);
}
QP.CPP
double* QuadEquLag::optimizer(double** H,int& m, int& n, double* c, double** A, double* b)
{
//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
double* res = new double[n];
inverseMatrix(H, n, n);
double** R = new double*[n];
for(int i = 0; i < n; i++)
{
R[i] = new double[m];
}
calculateR(H, m, n, A,R);
rTmp = R;
double** Q = new double*[n];
for(int i = 0; i < n; i++)
{
Q[i] = new double[n];
}
calculateQ(H, m, n, A,Q);
double* Qc = new double[n];
for(int i = 0; i < n; i++)
{
double tmp1 = 0;
double tmp2 = 0;
for(int k1 = 0; k1 < n; k1++)
{
tmp1 = tmp1 + Q[i][k1]*c[k1];
}
for(int k2 = 0; k2 < m; k2++)
{
tmp2 = tmp2 + R[i][k2]*b[k2];
}
res[i] = -tmp1 + tmp2;
}
for(int i = 0; i < n; i++)
{
delete[] R[i];
}
delete[] R;
for(int i = 0; i < n; i++)
{
delete[] Q[i];
}
delete[] Q;
return res;
}
void QuadEquLag::calculateQ(double** reverseH, int& m, int& n, double** A,double** Q)
{
//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
////================计算H逆乘A转置=====================
//double** rHtA = new double*[n]; //n*m维
//for(int i = 0; i < n; i++)
//{
// rHtA[i] = new double[m];
//}
//twoMatrixMultiply(reverseH,n,n,"N",A,m,n,"T",rHtA);
////=====================================================
////===============计算A乘H逆乘A转置=====================
//double** ArH = new double*[m]; //m*n维
//for(int i = 0; i < m; i++)
//{
// ArH[i] = new double[n];
//}
//twoMatrixMultiply(A,m,n,"N",reverseH,m,n,"N",ArH);
//double** ArHtA = new double*[m];
//for(int i = 0; i < m; i++)
//{
// ArHtA[i] = new double[m];
//}
//twoMatrixMultiply(ArH,m,n,"N",A,m,n,"T",ArHtA);
//inverseMatrix(ArHtA, m, m);
////释放内存空间
//for(int i = 0; i < m; i++)
//{
// delete[] ArH[i];
//}
//delete[] ArH;
//double** tmp1 = new double*[n];
//for(int i = 0; i < n; i++)
//{
// tmp1[i] = new double[m];
//}
//calculateR(reverseH, m, n, A,tmp1);
//=====================================================
//=====================计算A乘H逆======================
double** AtH = new double*[m];
for(int i = 0; i < m; i++)
{
AtH[i] = new double[n];
}
twoMatrixMultiply(A,m,n,"N", reverseH,n,n,"N",AtH);
//=====================================================
////释放内存
//for(int i = 0; i < n; i++)
//{
// delete[] rHtA[i];
//}
//delete[] rHtA;
//for(int i = 0; i < m; i++)
//{
// delete[] ArHtA[i];
//}
//delete[] ArHtA;
double** tmp2 = new double*[n];
for(int i = 0; i < n; i++)
{
tmp2[i] = new double[n];
}
twoMatrixMultiply(rTmp,n,m,"N", AtH,m,n,"N",tmp2);
//释放内存空间
/*for(int i =0; i < n; i++)
{
delete[] tmp1[i];
}
delete[] tmp1;*/
for(int i =0; i < m; i++)
{
delete[] AtH[i];
}
delete[] AtH;
//两个矩阵相减
for(int i = 0; i < n; i++)
{
for(int j = 0; j < n; j++)
{
Q[i][j] = reverseH[i][j] - tmp2[i][j];
}
}
//=================内存清理==================
for(int i =0; i < n; i++)
{
delete[] tmp2[i];
}
delete[] tmp2;
//===========================================
}
void QuadEquLag::calculateR(double** reverseH, int& m, int& n, double** A,double** R)
{
//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
//================计算H逆乘A转置=====================
double** rHtA = new double*[n]; //n*m维
for(int i = 0; i < n; i++)
{
rHtA[i] = new double[m];
}
twoMatrixMultiply(reverseH,n,n,"N",A,m,n,"T",rHtA);
//=====================================================
//===============计算A乘H逆乘A转置=====================
double** ArH = new double*[m]; //m*n维
for(int i = 0; i < m; i++)
{
ArH[i] = new double[n];
}
twoMatrixMultiply(A,m,n,"N",reverseH,n,n,"N",ArH);
double** ArHtA = new double*[m];
for(int i = 0; i < m; i++)
{
ArHtA[i] = new double[m];
}
twoMatrixMultiply(ArH,m,n,"N",A,m,n,"T",ArHtA);
inverseMatrix(ArHtA, m, m);
//释放内存空间
for(int i = 0; i < m; i++)
{
delete[] ArH[i];
}
delete[] ArH;
twoMatrixMultiply(rHtA,n,m,"N",ArHtA,m,m,"N",R);
//释放内存
for(int i = 0; i < n; i++)
{
delete[] rHtA[i];
}
delete[] rHtA;
for(int i = 0; i < m; i++)
{
delete[] ArHtA[i];
}
delete[] ArHtA;
}
void QuadEquLag::twoMatrixMultiply(double** a, int& aRow, int& aCol, char* aType, double** b, int& bRow, int& bCol, char* bType,double** out)
{
if(strcmp(aType,"N") == 0 && strcmp(bType,"N") == 0)//a*b
{
for(int i = 0; i < aRow; i++)
{
for(int j = 0; j < bCol; j++)
{
double tmp = 0;
for(int k = 0; k < aCol; k++)
{
tmp = tmp + a[i][k]*b[k][j];
}
out[i][j] = tmp;
}
}
}
else if(strcmp(aType,"N") == 0 && strcmp(bType,"T") == 0)//a*trans(b)
{
for(int i = 0; i < aRow; i++)
{
for(int j = 0; j < bRow; j++)
{
double tmp = 0.0;
for(int k = 0; k < bCol; k++)
{
tmp = tmp + a[i][k]*b[j][k];
}
out[i][j] = tmp;
}
}
}
else if(strcmp(aType,"T") == 0 && strcmp(bType,"N") == 0)//trans(a)*b
{
}
else if(strcmp(aType,"T") == 0 && strcmp(bType,"T") == 0)//trans(a)*trans(b)
{
}
}
void QuadEquLag::inverseMatrix(double** a, int& aRow, int& aCol)
{
int N = aRow;
double* A = new double[aRow * aCol];
int LDA = aCol;
int* IPIV = new int[N];
int LWORK = N;
double* WORK = new double[N];
int INFO ;
int k = 0;
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
A[k++] = a[j][i];
}
}
int m = N; int n = m;
dgetrf_(&m, &n, A, &m, IPIV, &INFO);
if(INFO != 0)
{
printf("dgetrf is wrong");
assert(0);
}
INFO = 1;
dgetri_(&N, A, &LDA, IPIV, WORK, &LWORK, &INFO);
if(INFO != 0)
{
printf("dgetri is wrong");
assert(0);
}
k = 0;
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
a[j][i] = A[k++];
}
}
delete[] A;
}
//double* QuadProg::optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0)
//{
// //H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
//
//
//
// //vector<int> J0 ;
// //
// ////1)计算有效集J0，构造矩阵Aa0，设k=0
// //for(int i = 0; i < m; i++)
// //{
// //
// // double tmp = 0.0;
// // for(int j = 0; j < n; j++)
// // {
// // tmp = tmp + A[i][j] * x0[j];
// // }
// // if(tmp == b[i])
// // {
// // J0.push_back(i);
// // }
// //}
// //int nn=0;
// //while(true)
// //{
// // //2)计算gk，并判断rank[trans(Aak) gk] 时候等于 rank(trans(Aak))
// // double* gk = new double[n];
// // for(int i = 0; i < n; i++)
// // {
// // double tmp = 0.0;
// // for(int j = 0; j < n; j++)
// // {
// // tmp = tmp + H[i][j] * x0[j];
// // }
// // gk[i] = c[i] + tmp;
// // }
// // double** Aakgk = new double*[n];
// // for(int i = 0; i < n; i++)
// // {
// // Aakgk[i] = new double[J0.size() + 1];
// // }
// //
// // for(unsigned int i = 0; i < J0.size(); i++)
// // {
// // for(int j = 0; j < n; j++)
// // {
// // Aakgk[j][i] = A[J0[i]][j];
// // }
// // }
// //
// // int last = J0.size();
// // for(int j = 0; j < n; j++)
// // {
// // Aakgk[j][last] = gk[j];
// // }
// //
// // int col = J0.size() + 1;
// // int rankAakgk = rank(Aakgk,n,col);
// // col--;
// // int rankAak = rank(Aakgk, n, col);
//
// // if(rankAakgk == rankAak)
// // {
// // //3)如果gk中的值全大于零，则迭代停止，找到最优值，gk相当于罚参数
// // solveLinearEquations(Aakgk, n, col, gk);//gk作为输入参数和输出参数，最终的方程组的解存储在此
// //
// // for(int i = 0; i < n; i++)
// // {
// // delete[] Aakgk[i];
// // }
// // delete[] Aakgk;
// //
// // /*if(greaterThanZero(gk,n) == true)
// // {
// // delete[] gk;
// // break;
// // }*/
// // if(greaterThanZero(gk,n) == false)
// // {
// // updateActiveSet(J0, gk );//更新有效集，将gk中的负值对应的有效集去掉，如果gk中全部为负值，去掉最小的负值所对应的有效集
// // /*delete[] gk;*/
// // }
// // }
// //
// // //4)求解dk
//
// // double** Aak = new double*[n];
// // int mTmp ;
// // double* bTmp = NULL;
// // if(J0.size() == 0)
// // {
// // for(int i = 0; i < n; i++)
// // {
// // Aak[i] = new double[1];
// // Aak[i][0] = 0.0;
// // }
// // mTmp = 1;
// // bTmp = new double[1];
// // for(unsigned int i = 0; i < 1; i++)
// // {
// // bTmp[i] = 0.0;
// // }
// // }
// // else
// // {
// // for(int i = 0; i < n; i++)
// // {
// // Aak[i] = new double[J0.size()];
// // }
// // for(unsigned int i = 0; i < J0.size(); i++)
// // {
// // for(int j = 0; j < n; j++)
// // {
// // Aak[j][i] = A[J0[i]][j];
// // }
// // }
// //
// // mTmp =(unsigned) J0.size();
// // bTmp = new double[J0.size()];
// // for(unsigned int i = 0; i < J0.size(); i++)
// // {
// // bTmp[i] = 0.0;
// // }
// // }
// // cout<<nn++<<endl;
//
// // double* dk = QuadEquLag::optimizer(H, mTmp, n, gk, Aak, bTmp);
// // if(equalZero(dk,n) == true && greaterThanZero(gk,n) == true)
// // {
// // delete[] gk;
// // break;
// // }
// // {
// // delete[] gk;
// // }
// // for(int i = 0; i < n; i++)
// // {
// // delete[] Aak[i];
// // }
// // delete[] Aak;
// // delete[] bTmp;
// // //5)求解afak，并设置x(k+1) = x(k) + afak*dk;
// // double* xk = x0;
// // int minIndex;
// // double afak = findAfak(J0, dk,A, b, xk,minIndex, m,n);
// // updataXk(xk,afak,dk,n);
//
// // //6)
// // if(afak < 1)
// // {
// // J0.push_back(minIndex);
// // }
// // sort(J0.begin(), J0.end());
// //}
//
// return x0;
//}
double QuadProg::optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0)
{
//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
vector<int> J0 ;
bool come5 = true;
//0)计算有效集J0，构造矩阵Aa0，设k=0
for(int i = 0; i < m; i++)
{
double tmp = 0.0;
for(int j = 0; j < n; j++)
{
tmp = tmp + A[i][j] * x0[j];
}
if(tmp == b[i])
{
J0.push_back(i);
}
}
int nn=0;
double* gk = new double[n];
int nm = 0;
while(true)
{
//1)计算gk，并判断rank[trans(Aak) gk] 时候等于 rank(trans(Aak))
cout << nm++ <<endl;
if(come5 == true)
{
come5 = false;
for(int i = 0; i < n; i++)
{
double tmp = 0.0;
for(int j = 0; j < n; j++)
{
tmp = tmp + H[i][j] * x0[j];
}
gk[i] = c[i] + tmp;
}
}
//2)求解线性方程组[Q A'; A 0]*[dk v]' = [-gk 0]
double* gkk = new double[n+J0.size()];
double** QA = new double*[(n+J0.size())];
for(unsigned int i = 0; i < n+J0.size(); i++)
{
QA[i] = new double[n+J0.size()];
}
for(int i = 0; i < n; i++)
{
for(int j = 0; j < n; j++)
{
QA[i][j] = H[i][j];
}
}
for(unsigned int i = 0; i < J0.size(); i++)
{
for(int j = 0; j < n; j++)
{
QA[i+n][j] = A[J0[i]][j];
}
}
for(unsigned int i = 0; i < J0.size(); i++)
{
for(int j = 0; j < n; j++)
{
QA[j][i+n] = A[J0[i]][j];
}
}
for(unsigned int i = 0; i < J0.size(); i++)
{
for(unsigned int j = 0; j < J0.size(); j++)
{
QA[i+n][j+n] = 0.0;
}
}
for(unsigned int i = 0; i < n+J0.size(); i++)
{
if(i < n)
gkk[i] = -gk[i];
else
gkk[i] = 0.0;
}
int row = n + J0.size();
int col = n + J0.size();
double* gkkk = new double[n+J0.size()];
for(unsigned int i = 0; i < n + J0.size(); i++)
{
gkkk[i] = gkk[i];
}
solveLinearEquations(QA, row, col, gkkk);
for(unsigned int i = 0; i < n+J0.size(); i++)
{
delete[] QA[i];
}
delete[] QA;
//3)如果dk=0,转到6）
if(equalZero(gkkk,n) != true)
{
//4)计算步长
int minIndex;
double afak =findAfak(J0, gkkk, A , b, x0, minIndex, m, n);
if(afak < 1)
{
J0.push_back(minIndex);
}
sort(J0.begin(), J0.end());
//5)更新xk
updataXk(x0,afak,gkkk,n);
come5 = true;
delete[] gkk;
delete[] gkkk;
continue;
}
else if(equalZero(gkkk,n) == true)
{
//6）
double* v = new double[J0.size()];
for(unsigned int i = 0; i < J0.size(); i++)
{
v[i] = gkkk[n+i];
}
int num = J0.size();
if(greaterThanZero(v,num) == true)
{
break;
}
else
{
updateActiveSet(J0, v);
}
delete[] v;
delete[] gkkk;
continue;
}
}
delete[] gk;
double obj = 0;
double* o = new double[n];
for(int i = 0; i < n; i++)
{
double tmp = 0.0;
for(int j = 0; j < n; j++)
{
tmp = tmp + x0[j]*H[i][j];
}
o[i] = tmp;
}
double cx = 0.0;
for(int i = 0; i < n; i++)
{
obj = obj + o[i] * x0[i] ;
cx = cx + c[i] * x0[i];
}
delete[] o;
return obj * 0.5 + cx;
;
}
bool QuadProg::equalZero(double* in, int& n)
{
for(int i = 0; i < n; i++)
{
if(in[i] > 0.0000000001)
return false;
}
return true;
}
void QuadProg::updataXk(double* xk,double& afak, double* dk, int& n)
{
for(int i = 0; i < n; i++)
{
xk[i] = xk[i] + afak * dk[i];
}
}
double QuadProg::findAfak(vector<int>& activeSet, double* dk, double** A, double* b, double* xk, int& minIndex, int& m, int& n)
{
double minVal = 0.0;
vector<double> minVals;
vector<int> minIndexs;
for(int i = 0; i < m; i++)
{
vector<int>::iterator iter = find(activeSet.begin(),activeSet.end(),i);
if(iter == activeSet.end()) //i不在有效集中
{
double aixk = 0.0;
double aidk = 0.0;
for(int j = 0; j < n; j++)
{
aidk = aidk + A[i][j] * dk[j];
aixk = aixk + A[i][j] * xk[j];
}
if(aidk > 0.0)
{
double rating = (-aixk + b[i])/(aidk);
minVals.push_back(rating);
minIndexs.push_back(i);
}
}
}
minVal = minVals[0];
minIndex = minIndexs[0];
for(unsigned int i = 1; i< minVals.size(); i++)
{
if(minVals[i] < minVal)
{
minVal = minVals[i];
minIndex = minIndexs[i];
}
}
if(1 < minVal)
return 1.0;
else
return minVal;
}
void QuadProg::updateActiveSet(vector<int>& activeSet, double* mu)
{
int m =(int) activeSet.size();
int minIndex = 0;
double minVal = mu[0];
for(int i = 1; i < m; i++)
{
if(mu[i] < minVal)
{
minVal = mu[i];
minIndex = i;
}
}
if(minVal < 0.0)
{
activeSet.erase(activeSet.begin() + minIndex);
}
}
bool QuadProg::greaterThanZero(double* in, int& n)
{
for(int i = 0; i < n; i++)
{
if(in[i] < 0.0)
return false;
}
return true;
}
int QuadProg::rank(double** a,int& row, int& col, double tol)
{
/*s = svd(A);rank = sum(s > tol);*/
//===========SVD=============
char jobu[]="A";
char jobvt[]="A";
int m = row;
int n = col;
double *b=new double[m*n];
for(int i=0;i<n;i++)
{
for(int j=0;j<m;j++)
{
b[i*m+j]=a[j][i];
}
}
int lda=m;
int minVal = min(m,n);
double *s=new double[minVal];
int ldu=m;
double *u=new double[ldu*m];
int ldvt=n;
double *vt=new double[ldvt*n];
int lwork=max(1,5*min(m,n));
double *work=new double[lwork];
int info;
dgesvd_(jobu, jobvt, &m, &n, b, &lda, s, u, &ldu, vt, &ldvt, work, &lwork,&info);
//===========================
delete[] b;
delete[] u;
delete[] vt;
int rank = 0;
for(int i = 0; i < minVal; i++)
{
if(s[i] > tol)
rank++;
}
delete[] s;
return rank;
}
void QuadProg::solveLinearEquations(double** a, int& row, int& col, double* b)
{
int N = row;
int NRHS = 1;
double* A = new double[(row) * (col)];
int index = 0;
for(int i = 0; i < row; i++)
{
for(int j = 0; j < col; j++)
{
A[index++] = a[j][i];
}
}
int LDA = row;
int* IPIV = new int[N];
double* B = b;
int LDB = max(1, N);
int INFO;
dgesv_(&N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO );
delete[] A;
delete[] IPIV;
}
测试主程序为：
#include "src\Optimizer\QP.h"
#include <fstream>
using namespace std;
int main(int argv, char** argc)
{
double** H = new double*[400];
double** A = new double*[400];
for(int i = 0; i < 400; i++)
{
A[i] = new double[400];
H[i] = new double[400];
}
ifstream readin;
readin.open("./data/H.txt");
if(readin.is_open())
{
for(int i = 0; i < 400; i++)
{
for(int j = 0; j < 400; j++)
{
double tmp;
readin >> tmp;
H[i][j] = tmp;
}
}
}
readin.close();
double* f = new double[400];
readin.open("./data/f.txt");
if(readin.is_open())
{
for(int i = 0; i < 400; i++)
{
double tmp;
readin >> tmp;
f[i] = tmp;
}
}
readin.close();
for(int i = 0; i < 400; i++)
{
for(int j = 0; j < 400; j++)
{
A[i][j] = 0.0;
}
A[i][i] = 1.0;
}
double* b = new double[400];
readin.open("./data/b.txt");
if(readin.is_open())
{
for(int i = 0; i < 400; i++)
{
double tmp;
readin >> tmp;
b[i] = tmp;
}
}
b[399] = 1000000000;
readin.close();
double* x0 = new double[400];
for(int i = 0; i < 400; i++)
{
x0[i] = 0.5;
}
int m = 400;
int n = 400;
QuadProg* qu = new QuadProg();
double obj = qu->optimizer( H,m,n, f, A, b, x0);
return 1;
}
本程序中用到了Lapack库和测试数据（Lapack是进行矩阵计算的库，Matlab中也采用了它），可以到此处下载，漏掉的Lapack.h文件见这里：点击打开链接