• wqz11234
    了解作者
  • matlab
    开发工具
  • 5KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 4
    下载次数
  • 2020-09-18 12:02
    上传日期
大规模SDP求解器,基于SDPT3实现文档,全部MATLAB实现,没有底层c语言库
SDPSlover.rar
  • SDPSlover
  • mat.m
    232B
  • svec.m
    611B
  • SDPSlover.m
    7.7KB
  • svec_1block.m
    363B
  • H.m
    284B
  • block.m
    350B
  • smat_1block.m
    419B
  • vblock.m
    239B
  • skmult.m
    483B
  • smat.m
    533B
内容介绍
% SDP Solver clc clear tic; %SDPA格式 data=[ 2 ... 2 ... 2 2 ... 10.0 20.0... 0 1 1 1 1.0... 0 1 2 2 2.0... 0 2 1 1 3.0... 0 2 2 2 4.0... 1 1 1 1 1.0... 1 1 2 2 1.0... 2 1 2 2 1.0... 2 2 1 1 5.0... 2 2 1 2 2.0... 2 2 2 2 6.0]; data = data'; global n ns nis n2 n2s n2is nt nts ntis nblocks; % m = 约束矩阵数量 m = data(1); % nblocks = 约束矩阵的块对角线结构中的块数 nblocks = data(2); % 1D 索引 % n's = 每个nxn块的1D大小(n)(0是前缀) ns = data(3:3+nblocks-1); ns = [ 0; abs(ns)]; % n = sum( n's ) n = sum(ns); % n_i's = 每个块的累积一维指数 nis = zeros(size(ns)); for i = 1:length(ns) nis(i) = sum(ns(1:i)); end % 2D 索引 % n^2's = 每个nxn块的2D大小(n*n)(0前缀) n2s = ns.^2; % n2 = sum( n^2's ) n2 = sum(n2s); % n^2_i's = 每个块的累积2D索引 n2is = zeros(size(ns)); for i = 1:length(ns) n2is(i) = sum(n2s(1:i)); end % svec 索引 % nt's = (n*(n+1)/2) svec的三角索引(nxn块) (0前缀) nts = ns.*(ns+1)/2; % nt_i's = svec的累积三角索引(nxn块) ntis = zeros(size(nts)); for i = 1:length(nts) ntis(i) = sum(nts(1:i)); end % nt = sum(nt's) = sum( n_i*(n_i+1)/2 ) nt = sum(nts); % 从文件中读取数据并分配空间 b = -data(3+nblocks:3+nblocks+m-1); data = data(3+nblocks+m:end); c = zeros(n2, 1); A = zeros(m, n2); % 读取约束矩阵 for idx = 1:5:size(data,1)-1 matno = data(idx); %第几个约束矩阵 blkno = data(idx+1); %第几个块 i = data(idx+2); %行 j = data(idx+3); %列 entry = data(idx+4);%值 id1 = (j-1)*ns(blkno+1) + i + n2is(blkno); %行索引 id2 = (i-1)*ns(blkno+1) + j + n2is(blkno); %列索引 if matno == 0 c(id1) = -entry; c(id2) = -entry; else A(matno, id1) = -entry; A(matno, id2) = -entry; end end % 将A和c转换为svec版本 c = mat(c); c = svec(c); A_temp = A; A = zeros(m, nt); for i = 1:m A(i,:) = svec(mat(A_temp(i,:))); end % 初始X,Z选取,基于SDPT3文档 % 设ksi ksi = zeros(nblocks,1); temp = zeros(m, 1); for j = 1:nblocks for k = 1:m %求svec A中某个块的所有元素的范数,再用1+b/1+normA块, 获得缩放系数ksi temp(k) = (1 + abs(b(k))) / (1 + norm(A(k, vblock(j)))); end ksi(j) = ns(j+1) * max(temp); end % 设缩放系数eta eta = zeros(nblocks, 1); temp = zeros(m,1); for j = 1:nblocks for k = 1:m temp(k) = norm(A(k, vblock(j)));%求svec A中某个块的所有元素的范数 end eta(j) = (1/sqrt(ns(j+1))) * (1 + max( max(temp), norm(c(vblock(j))))); end % 设置初始迭代使用的X和Z X = zeros(n); for j = 1:nblocks X(block(j)) = eye(ns(j+1))*ksi(j); end x = svec(X); y = zeros(m,1); Z = zeros(n); for j = 1:nblocks Z(block(j)) = eye(ns(j+1))*eta(j); end z = svec(Z); % 计算每个块的normA,保存每个块的sqrt normsA = zeros(nblocks,1); for j = 1:nblocks normsA(j) = norm(A(:, vblock(j)), 'fro');%矩阵的各个元素绝对值的最大值 end normsA = max(1, sqrt(normsA)); % 计算normb normb = max(1, norm(b)); % 计算每个块 normc , 保留max normc = 1; for j = 1:nblocks normc = max(normc, norm(c(vblock(j)))); end % 缩放约束矩阵与系数向量 for j = 1:nblocks A(:, vblock(j)) = A(:, vblock(j)) / normsA(j); c(vblock(j)) = c(vblock(j)) / (normc*normsA(j)); end b = b / normb; % 缩放初始迭代x和z for j = 1:nblocks x(vblock(j)) = x(vblock(j)) * normsA(j); z(vblock(j)) = z(vblock(j)) / (normc*normsA(j)); end % 设置初始值 rp = b - A*x; Rd = c - z -A'*y; relgap = dot(x,z) / (1 + max(abs(dot(c,x)), abs(dot(b,y)))); %停止条件1:相对对偶间隙relgap pinfeas = norm(rp) / (1 + norm(b)); dinfeas = norm(Rd) / (1 + norm(c)); phi = max(pinfeas, dinfeas); %停止条件2:不可行系数phi soln_relgap = 1e-6; %停止准则 soln_phi = 1e-6; %停止准则 gamma = 0.9; iter = 0; % 更新中间解直到收敛 while(relgap > soln_relgap || phi > soln_phi) % 得到X和Z的Cholesky分解 X = smat(x); Q = chol(X); Z = smat(z); P = chol(Z); % 计算残差mu,relgap,phi %(4) rp = b - A*x; Rd = c - z -A'*y; %svec Rd mu = dot(x,z) / n; %mu relgap = dot(x,z) / (1 + max(abs(dot(c,x)), abs(dot(b,y)))); phi = max(norm(rp)/(1+norm(b)), norm(Rd)/(1+norm(c))); Rc = -svec( H(X*Z, P) ); % 计算预测方向(伪方向 sigma =0) (10-16) HKM方向 得到dx,dz的svec形式 M = A*skmult(X, Z\eye(n), A'); %NT (43) h = rp + A*(skmult(X, Z\eye(n), Rd) - skmult(P\eye(n), P\eye(n), Rc)); dy = M\h; dz = Rd - A'*dy; %dx = -x - skmult(X, Z\eye(n), dz); dx = skmult(P\eye(n), P\eye(n), Rc) - skmult(X, Z\eye(n), dz); % 计算 alpha (原始步长)(33) alphas = zeros(1,nblocks); for j = 1:nblocks XDX = (Q(block(j))') \ smat_1block(dx(vblock(j))) / Q(block(j));%(32) eigs_a = sort(eig(XDX)); %排序 lambda_a = eigs_a(1); %取最小的特征值 if lambda_a < 0 alphas(j) = -1/lambda_a; else alphas(j) = Inf; end end alpha = min(1, gamma*min(alphas)); % 计算 beta (对偶步长)(34) betas = zeros(1,nblocks); for j = 1:nblocks ZDZ = (P(block(j))') \ smat_1block(dz(vblock(j))) / P(block(j)); eigs_b = sort(eig(ZDZ)); lambda_b = eigs_b(1); if lambda_b < 0 betas(j) = -1/lambda_b; else betas(j) = Inf; end end beta = min(1, gamma*min(betas)); % 计算指数e (sigma的指数) if mu > 1e-6 if min(alpha, beta) < 1 / sqrt(3) e = 1; else e = max(1, 3*min(alpha, beta)^2); end else e = 1; end % 计算 sigma (中心参数) if dot(x+alpha*dx, z+beta*dz) < 0 sigma = 0.8; else frac = dot(x+alpha*dx, z+beta*dz)/dot(x,z); sigma = min(1, frac^e); end % 计算 正确 的预测方向 (10-16) HKM方向 得到dx,dz的svec形式 Rc = svec(sigma*mu*eye(n)-H(X*Z, P)-H(smat(dx)*smat(dz), P)); h = rp + A*(skmult(X, Z\eye(n), Rd) - skmult(P\eye(n), P\eye(n), Rc)); dy = M\h; dz = Rd - A'*dy; dx = skmult(P\eye(n), P\eye(n), Rc) - skmult(X, Z\eye(n), dz); gamma = 0.9 + 0.09*min(alpha, beta); % 计算正确 alpha (原始步长)(33) for j = 1:nblocks XDX = (Q(block(j))') \ smat_1block(dx(vblock(j))) / Q(block(j)); eigs_a = sort(eig(XDX)); lambda_a = eigs_a(1); if lambda_a < 0 alphas(j) = -1/lambda_a; else alphas(j) = Inf; end end alpha = min(1, gamma*min(alphas)); % 计算正确 beta (对偶步长)(34) for j = 1:nblocks ZDZ = (P(block(j))') \ smat_1block(dz(vblock(j))) / P(block(j)); eigs_b = sort(eig(ZDZ)); lambda_b = eigs_b(1); if lambda_b < 0 betas(j) = -1/lambda_b; else betas(j) = Inf; end end beta = min(1, gamma*min(betas)); % 更新步骤 x = x + alpha*dx; y = y + beta*dy; z = z + beta*dz; gamma = 0.9 + 0.09*min(alpha, beta); % 打印迭代结果 iter = iter + 1; primal = dot(c,x); %原始结果 dual = dot(b,y); %对偶结果 pinfeas = norm(rp) / (1 + norm(b)); %原始可行性 dinfeas = norm(Rd) / (1 + norm(c)); %对偶可行性 fprintf('%d\t%.3f\t%.3f\t%.3e\t%.3e %.3f\t%.2e %.2e %.2e %.2e %.2e %.2e %.2f %2e %2e\n', ... iter, alpha, beta, sigma, mu, e, pinfeas, dinfeas, relgap, primal, dual, phi, gamma, norm(rp), norm(Rd)); end disp(' '); disp('计算完毕'); % 反缩放约束矩阵 for j = 1:nblocks A(:, vblock(j)) = A(:, vblock(j)) * normsA(j); c(vblock(j)) = c(vblock(j)) * (normc*normsA(j)); end b = b * normb; % 反缩放结果 for j = 1:nblocks x(vblock(j)) = x(vblock(j)) * normb / normsA(j); z(vblock(j)) = z(vblock(j)) * (normc*normsA(j)); end y = y * normc; X = smat(x); Z = smat(z); % 打印结果 disp(['原始结果: ', num2str(dot(-c,x))]); disp(['对偶结果: ', num2str(dot(-b,y))]); time = toc; disp(['时间: ', num2str(time), ' s']);
评论
    相关推荐
    • matlabcnhelp.rar
      matlab中文帮助很难找的,快速下载
    • MobilePolice.rar
      移动警察,车牌识别,车牌定位系统源代码,已经运用在移动车载稽查系统中。
    • SVM(matlab).rar
      支持向量机(SVM)实现的分类算法源码[matlab]
    • svm.zip
      用MATLAB编写的svm源程序,可以实现支持向量机,用于特征分类或提取
    • Classification-MatLab-Toolbox.rar
      模式识别matlab工具箱,包括SVM,ICA,PCA,NN等等模式识别算法,很有参考价值
    • VC++人脸定位实例.rar
      一个经典的人脸识别算法实例,提供人脸五官定位具体算法及两种实现流程.
    • QPSK_Simulink.rar
      QPSK的Matlab/Simulink的调制解调仿真系统,给出接收信号眼图及系统仿真误码率,包含载波恢复,匹配滤波,定时恢复等重要模块,帮助理解QPSK的系统
    • LPRBPDemo2009KV.rar
      车牌识别,神经网络算法,识别率高达95%,识别时间低于80ms。
    • MODULATION.RAR
      这个源程序代码包提供了通信系统中BPSK,QPSK,OQPSK,MSK,MSK2,GMSK,QAM,QAM16等调制解调方式 用matlab的实现,以及它们在AWGN和Rayleigh信道下的通信系统实现及误码率性能
    • algorithms.rar
      十大算法论文,包括遗传算法,模拟退火,蒙特卡罗法等等,对于初学者很有帮助!!