Skip to content

Instantly share code, notes, and snippets.

@chinasaur

chinasaur/Handcoded

Created September 13, 2013 04:45
Show Gist options
  • Save chinasaur/6546857 to your computer and use it in GitHub Desktop.
Save chinasaur/6546857 to your computer and use it in GitHub Desktop.
Download ZIP
Matlab profiling comparison for SOCP program stuffing
Raw
Handcoded
time calls line
1 function [ecos_c, ecos_G, ecos_h, dims, ecos_A, ecos_b] = cbp2ecos_2norm(dict, data, lambda, noisesigma, radii, theta)
2 % RADII, THETA must be column vectors
3
1000 4 ndict = size(dict, 2) / 3;
1000 5 ldata = length(data);
6
7
8 % Cone constraints in conventional order (ecos_G in particular will have to
9 % be arranged to match these conventions)
1000 10 dims.f = ndict;
1000 11 dims.l = ndict*2;
1000 12 dims.q = [ldata+1; 3.*ones(ndict,1)];
13
14
15 % Calculate locations for variables
1000 16 slack1 = 1;
1000 17 slack2 = ndict*2;
1000 18 quadobj = slack2+1;
1000 19 slack3 = quadobj+1;
< 0.01 1000 20 slack4 = quadobj+ldata;
1000 21 quadconstr1 = slack4+1;
1000 22 quadconstr2 = slack4+ndict*3;
23
24
1000 25 ecos_c = data;
1000 26 ecos_c = [ecos_c; zeros(ndict*2,1)];
< 0.01 1000 27 ecos_h = zeros(quadconstr2,1);
< 0.01 1000 28 ecos_h(quadobj) = 1 / sqrt(2) / noisesigma;
29
30
31 % ecos_G...
32
33 % Objective norm constraint slacks
1000 34 Oncs = -speye(ldata);
35
36 % Linear constraint slacks
0.04 1000 37 Lcs = -speye(ndict*2);
38
39 % Radii constraint dictionary equalities
0.24 1000 40 Rcde = reshape([sparse(ldata,ndict); -dict(:,2:3:end); -dict(:,3:3:end)], ldata, ndict*3)';
41
42 % Radii constraint eye
< 0.01 1000 43 Rceye1 = reshape([-speye(ndict); sparse(ndict*2,ndict)], ndict, ndict*3)';
0.04 1000 44 Rceye2 = reshape([ sparse(ndict,ndict); speye(ndict); sparse(ndict,ndict)], ndict, ndict*3)';
< 0.01 1000 45 Rceye = [Rceye1 Rceye2];
46
0.04 1000 47 ecos_G1 = [sparse(ndict*2,ldata); sparse(1,ldata); Oncs; Rcde ];
0.02 1000 48 ecos_G2 = [Lcs ; sparse(1,ndict*2); sparse(ldata,ndict*2); Rceye];
0.05 1000 49 ecos_G = [ecos_G1 ecos_G2];
50
51
52 % Free variables (c_coeffs)
0.11 1000 53 F = sparse(-dict(:,1:3:end)');
54
55 % Radii constraints
< 0.01 1000 56 Rc = sparse(1:ndict, 1:ndict, radii);
57
58 % Linear constraints
0.02 1000 59 Lc = sparse(1:ndict, 1:ndict, -radii.*cos(theta));
60
0.02 1000 61 ecos_A = [F Rc Lc];
< 0.01 1000 62 ecos_b = lambda;
Raw
QCML
time calls line
1 function [data] = prob_to_socp(params)
2 % PROB2SOCP: maps PARAMS into a struct of SOCP matrices
3 % Where input struct PARAMS has the following fields:
4 % 'dictu' has shape Matrix(83,56)
5 % 'noise' has shape Scalar()
6 % 'rctheta' has shape Matrix(56,56)
7 % 'radii' has shape Matrix(56,56)
8 % 'dictv' has shape Matrix(83,56)
9 % 'dictc' has shape Matrix(83,56)
10 % 'data' has shape Vector(83)
11 % 'lambda' has shape Vector(56)
12
1000 13 p = 0; m = 364; n = 225;
1000 14 c = zeros(n,1);
1000 15 h = zeros(m,1);
1000 16 b = zeros(p,1);
< 0.01 1000 17 Gi = []; Gj = []; Gv = [];
1000 18 Ai = []; Aj = []; Av = [];
19
0.02 1000 20 dims.l = 112;
< 0.01 1000 21 dims.q = [84; 3*ones(56,1)];
1000 22 dims.s = [];
23
24 % stuffing the objective vector
1000 25 c(1:1) = params.noise;
1000 26 c(58:113) = params.lambda;
27
28 % for the constraint rctheta*c + -1*u <= 0
0.02 1000 29 Gi = [Gi; mod(find(params.rctheta)-1,size(params.rctheta,1))];
< 0.01 1000 30 Gj = [Gj; 57 + floor((find(params.rctheta)-1)/size(params.rctheta,1))];
1000 31 Gv = [Gv; nonzeros(params.rctheta)];
< 0.01 1000 32 Gi = [Gi; (0:55)'];
1000 33 Gj = [Gj; (113:168)'];
0.02 1000 34 Gv = [Gv; -1*ones(56,1)];
35
36 % for the constraint _t1 + -1*radii*c <= 0
1000 37 Gi = [Gi; (56:111)'];
1000 38 Gj = [Gj; (1:56)'];
1000 39 Gv = [Gv; 1*ones(56,1)];
0.03 1000 40 Gi = [Gi; 56 + mod(find(params.radii)-1,size(params.radii,1))];
0.03 1000 41 Gj = [Gj; 57 + floor((find(params.radii)-1)/size(params.radii,1))];
0.03 1000 42 Gv = [Gv; nonzeros(-params.radii)];
43
44 % for the SOC constraint norm([data + -1*dictc*c + -1*dictv*v + -1*dictu*u]) <= _t0
1000 45 h(114:196) = params.data;
0.12 1000 46 Gi = [Gi; 113 + mod(find(params.dictc)-1,size(params.dictc,1))];
0.08 1000 47 Gj = [Gj; 57 + floor((find(params.dictc)-1)/size(params.dictc,1))];
0.02 1000 48 Gv = [Gv; nonzeros(params.dictc)];
0.07 1000 49 Gi = [Gi; 113 + mod(find(params.dictu)-1,size(params.dictu,1))];
0.04 1000 50 Gj = [Gj; 113 + floor((find(params.dictu)-1)/size(params.dictu,1))];
0.02 1000 51 Gv = [Gv; nonzeros(params.dictu)];
0.08 1000 52 Gi = [Gi; 113 + mod(find(params.dictv)-1,size(params.dictv,1))];
0.08 1000 53 Gj = [Gj; 169 + floor((find(params.dictv)-1)/size(params.dictv,1))];
0.03 1000 54 Gv = [Gv; nonzeros(params.dictv)];
1000 55 Gi = [Gi; 112];
< 0.01 1000 56 Gj = [Gj; 0];
1000 57 Gv = [Gv; -1];
58
59 % for the SOC product constraint norm(u, v) <= _t1
1000 60 Gi = [Gi; (198:3:365)'];
< 0.01 1000 61 Gj = [Gj; (169:224)'];
0.03 1000 62 Gv = [Gv; -1*ones(56,1)];
1000 63 Gi = [Gi; (197:3:364)'];
1000 64 Gj = [Gj; (113:168)'];
< 0.01 1000 65 Gv = [Gv; -1*ones(56,1)];
0.02 1000 66 Gi = [Gi; (196:3:363)'];
< 0.01 1000 67 Gj = [Gj; (1:56)'];
< 0.01 1000 68 Gv = [Gv; -1*ones(56,1)];
69
70 % Convert from sparse triplet to column compressed format.
71 % Also convert from 0 indexed to 1 indexed.
1000 72 A = sparse(Ai+1, Aj+1, Av, p, n);
1.17 1000 73 G = sparse(Gi+1, Gj+1, Gv, m, n);
74
75 % Build output
0.04 1000 76 data = struct('c', c, 'b', b, 'h', h, 'G', G, 'A', A, 'dims', dims);
< 0.01 1000 77 end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment