totsu/problem/
sdp.rs

1use num_traits::{Float, Zero, One};
2use totsu_core::solver::{Solver, SliceLike, Operator, Cone};
3use totsu_core::{LinAlgEx, MatType, MatOp, ConePSD, ConeZero, splitm, splitm_mut};
4use crate::MatBuild;
5
6//
7
8pub struct ProbSDPOpC<'a, L: LinAlgEx>
9{
10    vec_c: MatOp<'a, L>,
11}
12
13impl<'a, L: LinAlgEx> Operator<L> for ProbSDPOpC<'a, L>
14{
15    fn size(&self) -> (usize, usize)
16    {
17        let (n, one) = self.vec_c.size();
18        assert_eq!(one, 1);
19
20        (n, 1)
21    }
22
23    fn op(&self, alpha: L::F, x: &L::Sl, beta: L::F, y: &mut L::Sl)
24    {
25        // y = a*vec_c*x + b*y;
26        self.vec_c.op(alpha, x, beta, y);
27    }
28
29    fn trans_op(&self, alpha: L::F, x: &L::Sl, beta: L::F, y: &mut L::Sl)
30    {
31        // y = a*vec_c^T*x + b*y;
32        self.vec_c.trans_op(alpha, x, beta, y);
33    }
34
35    fn absadd_cols(&self, tau: &mut L::Sl)
36    {
37        self.vec_c.absadd_cols(tau);
38    }
39
40    fn absadd_rows(&self, sigma: &mut L::Sl)
41    {
42        self.vec_c.absadd_rows(sigma);
43    }
44}
45
46//
47
48pub struct ProbSDPOpA<'a, L: LinAlgEx>
49{
50    symmat_f: MatOp<'a, L>,
51    mat_a: MatOp<'a, L>,
52}
53
54impl<'a, L: LinAlgEx> ProbSDPOpA<'a, L>
55{
56    fn dim(&self) -> (usize, usize, usize)
57    {
58        let (sk, n) = self.symmat_f.size();
59        let (p, n_) = self.mat_a.size();
60        assert_eq!(n, n_);
61
62        (n, sk, p)
63    }
64}
65
66impl<'a, L: LinAlgEx> Operator<L> for ProbSDPOpA<'a, L>
67{
68    fn size(&self) -> (usize, usize)
69    {
70        let (n, sk, p) = self.dim();
71
72        (sk + p, n)
73    }
74
75    fn op(&self, alpha: L::F, x: &L::Sl, beta: L::F, y: &mut L::Sl)
76    {
77        let (_n, sk, p) = self.dim();
78
79        splitm_mut!(y, (y_sk; sk), (y_p; p));
80
81        // y_sk = a*symmat_f*x + b*y_sk
82        self.symmat_f.op(alpha, x, beta, &mut y_sk);
83
84        // y_p = a*mat_a*x + b*y_p
85        self.mat_a.op(alpha, x, beta, &mut y_p);
86    }
87
88    fn trans_op(&self, alpha: L::F, x: &L::Sl, beta: L::F, y: &mut L::Sl)
89    {
90        let (_n, sk, p) = self.dim();
91
92        splitm!(x, (x_sk; sk), (x_p; p));
93
94        // y = a*symmat_f^T*x_sk + a*mat_a^T*x_p + b*y
95        self.symmat_f.trans_op(alpha, &x_sk, beta, y);
96        self.mat_a.trans_op(alpha, &x_p, L::F::one(), y);
97    }
98
99    fn absadd_cols(&self, tau: &mut L::Sl)
100    {
101        self.symmat_f.absadd_cols(tau);
102        self.mat_a.absadd_cols(tau);
103    }
104
105    fn absadd_rows(&self, sigma: &mut L::Sl)
106    {
107        let (_n, sk, p) = self.dim();
108
109        splitm_mut!(sigma, (sigma_sk; sk), (sigma_p; p));
110
111        self.symmat_f.absadd_rows(&mut sigma_sk);
112        self.mat_a.absadd_rows(&mut sigma_p);
113    }
114}
115
116//
117
118pub struct ProbSDPOpB<'a, L: LinAlgEx>
119{
120    symvec_f_n: MatOp<'a, L>,
121    vec_b: MatOp<'a, L>,
122}
123
124impl<'a, L: LinAlgEx> ProbSDPOpB<'a, L>
125{
126    fn dim(&self) -> (usize, usize, usize)
127    {
128        let (sk, n) = self.symvec_f_n.size();
129        let (p, one) = self.vec_b.size();
130        assert_eq!(one, 1);
131
132        (n, sk, p)
133    }
134}
135
136impl<'a, L: LinAlgEx> Operator<L> for ProbSDPOpB<'a, L>
137{
138    fn size(&self) -> (usize, usize)
139    {
140        let (_n, sk, p) = self.dim();
141
142        (sk + p, 1)
143    }
144
145    fn op(&self, alpha: L::F, x: &L::Sl, beta: L::F, y: &mut L::Sl)
146    {
147        let (_n, sk, p) = self.dim();
148
149        splitm_mut!(y, (y_sk; sk), (y_p; p));
150
151        // y_sk = a*-symmat_f*x + b*y_sk
152        self.symvec_f_n.op(-alpha, x, beta, &mut y_sk);
153
154        // y_p = a*vec_b*x + b*y_p
155        self.vec_b.op(alpha, x, beta, &mut y_p);
156    }
157
158    fn trans_op(&self, alpha: L::F, x: &L::Sl, beta: L::F, y: &mut L::Sl)
159    {
160        let (_n, sk, p) = self.dim();
161
162        splitm!(x, (x_sk; sk), (x_p; p));
163
164        // y = a*-symvec_f_n^T*x_sk + a*vec_b^T*x_p + b*y
165        self.symvec_f_n.trans_op(-alpha, &x_sk, beta, y);
166        self.vec_b.trans_op(alpha, &x_p, L::F::one(), y);
167    }
168
169    fn absadd_cols(&self, tau: &mut L::Sl)
170    {
171        self.symvec_f_n.absadd_cols(tau);
172        self.vec_b.absadd_cols(tau);
173    }
174
175    fn absadd_rows(&self, sigma: &mut L::Sl)
176    {
177        let (_n, sk, p) = self.dim();
178
179        splitm_mut!(sigma, (sigma_sk; sk), (sigma_p; p));
180
181        self.symvec_f_n.absadd_rows(&mut sigma_sk);
182        self.vec_b.absadd_rows(&mut sigma_p);
183    }
184}
185
186//
187
188pub struct ProbSDPCone<'a, L: LinAlgEx>
189{
190    sk: usize,
191    p: usize,
192    cone_psd: ConePSD<'a, L>,
193    cone_zero: ConeZero<L>,
194}
195
196impl<'a, L: LinAlgEx> Cone<L> for ProbSDPCone<'a, L>
197{
198    fn proj(&mut self, dual_cone: bool, x: &mut L::Sl) -> Result<(), ()>
199    {
200        let (sk, p) = (self.sk, self.p);
201
202        splitm_mut!(x, (x_sk; sk), (x_p; p));
203
204        self.cone_psd.proj(dual_cone, &mut x_sk)?;
205        self.cone_zero.proj(dual_cone, &mut x_p)?;
206        Ok(())
207    }
208
209    fn product_group<G: Fn(&mut L::Sl) + Copy>(&self, dp_tau: &mut L::Sl, group: G)
210    {
211        let (sk, p) = (self.sk, self.p);
212
213        splitm_mut!(dp_tau, (t_sk; sk), (t_p; p));
214
215        self.cone_psd.product_group(&mut t_sk, group);
216        self.cone_zero.product_group(&mut t_p, group);
217    }
218    
219}
220
221//
222
223/// Semidefinite program
224/// 
225/// <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
226/// <script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-svg.js"></script>
227/// 
228/// The problem is
229/// \\[
230/// \begin{array}{ll}
231/// {\rm minimize} & c^Tx \\\\
232/// {\rm subject \ to} & \sum_{i=0}^{n - 1} x_i F_i + F_n \preceq 0 \\\\
233/// & A x = b,
234/// \end{array}
235/// \\]
236/// where
237/// - variables \\( x \in \mathbb{R}^n \\)
238/// - \\( c \in \mathbb{R}^n \\)
239/// - \\( F_j \in \mathcal{S}^k \\) for \\( j = 0, \ldots, n \\)
240/// - \\( A \in \mathbb{R}^{p \times n},\ b \in \mathbb{R}^p \\).
241/// 
242/// This is already a conic problem and can be reformulated as follows:
243/// \\[
244/// \begin{array}{ll}
245/// {\rm minimize} & c^Tx \\\\
246/// {\rm subject \ to} &
247///   \left[ \begin{array}{ccc}
248///   {\rm vec}(F_0) & \cdots & {\rm vec}(F_{n - 1}) \\\\
249///   & A &
250///   \end{array} \right]
251///   x + s =
252///   \left[ \begin{array}{c}
253///   -{\rm vec}(F_n) \\\\ b
254///   \end{array} \right] \\\\
255/// & s \in {\rm vec}(\mathcal{S}_+^k) \times \lbrace 0 \rbrace^p.
256/// \end{array}
257/// \\]
258/// 
259/// \\( {\rm vec}(X) = (X_{11}\ \sqrt2 X_{12}\ X_{22}\ \sqrt2 X_{13}\ \sqrt2 X_{23}\ X_{33}\ \cdots)^T \\)
260/// which extracts and scales the upper-triangular part of a symmetric matrix X in column-wise.
261/// [`ConePSD`] is used for \\( {\rm vec}(\mathcal{S}_+^k) \\).
262pub struct ProbSDP<L: LinAlgEx>
263{
264    vec_c: MatBuild<L>,
265    mat_a: MatBuild<L>,
266    vec_b: MatBuild<L>,
267
268    symmat_f: MatBuild<L>,
269    symvec_f_n: MatBuild<L>,
270
271    eps_zero: L::F,
272    w_cone_psd: Vec<L::F>,
273    w_solver: Vec<L::F>,
274}
275
276impl<L: LinAlgEx> ProbSDP<L>
277{
278    /// Creates a SDP with given data.
279    /// 
280    /// Returns a [`ProbSDP`] instance.
281    /// * `vec_c` is \\(c\\).
282    /// * `syms_f` is \\(F_0, \\ldots, F_n\\) each of which shall belong to [`MatType::SymPack`].
283    /// * `mat_a` is \\(A\\).
284    /// * `vec_b` is \\(b\\).
285    /// * `eps_zero` should be the same value as [`totsu_core::solver::SolverParam::eps_zero`].
286    pub fn new(
287        vec_c: MatBuild<L>,
288        mut syms_f: Vec<MatBuild<L>>,
289        mat_a: MatBuild<L>, vec_b: MatBuild<L>,
290        eps_zero: L::F) -> Self
291    {
292        let n = vec_c.size().0;
293        let p = vec_b.size().0;
294
295        assert_eq!(vec_c.size(), (n, 1));
296        assert_eq!(syms_f.len(), n + 1);
297        let k = syms_f[0].size().0;
298        for sym_f in syms_f.iter() {
299            assert!(sym_f.is_sympack());
300            assert_eq!(sym_f.size(), (k, k));
301        }
302        assert_eq!(mat_a.size(), (p, n));
303        assert_eq!(vec_b.size(), (p, 1));
304
305        let f1 = L::F::one();
306        let f2 = f1 + f1;
307        let fsqrt2 = f2.sqrt();
308    
309        for sym_f in syms_f.iter_mut() {
310            sym_f.set_scale_nondiag(fsqrt2);
311            sym_f.set_reshape_colvec();
312        }
313
314        let symvec_f_n = syms_f.pop().unwrap();
315        let sk = symvec_f_n.size().0;
316
317        let symmat_f = MatBuild::new(MatType::General(sk, n))
318                       .by_fn(|r, c| { syms_f[c][(r, 0)] });
319
320        ProbSDP {
321            vec_c,
322            mat_a,
323            vec_b,
324            symmat_f,
325            symvec_f_n,
326            eps_zero,
327            w_cone_psd: Vec::new(),
328            w_solver: Vec::new(),
329        }
330    }
331
332    /// Generates the problem data structures to be fed to [`Solver::solve`].
333    /// 
334    /// Returns a tuple of operators, a cone and a work slice.
335    pub fn problem(&mut self) -> (ProbSDPOpC<L>, ProbSDPOpA<L>, ProbSDPOpB<L>, ProbSDPCone<'_, L>, &mut[L::F])
336    {
337        let p = self.vec_b.size().0;
338        let sk = self.symvec_f_n.size().0;
339
340        let f0 = L::F::zero();
341
342        let op_c = ProbSDPOpC {
343            vec_c: self.vec_c.as_op(),
344        };
345        let op_a = ProbSDPOpA {
346            symmat_f: self.symmat_f.as_op(),
347            mat_a: self.mat_a.as_op(),
348        };
349        let op_b = ProbSDPOpB {
350            symvec_f_n: self.symvec_f_n.as_op(),
351            vec_b: self.vec_b.as_op(),
352        };
353
354        self.w_cone_psd.resize(ConePSD::<L>::query_worklen(sk), f0);
355        let cone = ProbSDPCone {
356            sk, p,
357            cone_psd: ConePSD::new(self.w_cone_psd.as_mut(), self.eps_zero),
358            cone_zero: ConeZero::new(),
359        };
360
361        self.w_solver.resize(Solver::<L>::query_worklen(op_a.size()), f0);
362
363        (op_c, op_a, op_b, cone, self.w_solver.as_mut())
364    }
365}
totsu/problem/sdp.rs

totsu/problem/
sdp.rs
