ellalgo_rs/oracles/ldlt_mgr.rs
1use ndarray::{s, Array1, Array2};
2// use ndarray_linalg::Lapack;
3// use std::cmp::min;
4
5/// The `LDLTMgr` struct is a manager for LDLTMgr factorization in Rust.
6///
7/// `LDLTMgr` is a class that performs the LDLTMgr factorization for a given
8/// symmetric matrix. The LDLTMgr factorization decomposes a symmetric matrix A into
9/// the product of a lower triangular matrix L, a diagonal matrix D, and the
10/// transpose of L. This factorization is useful for solving linear systems and
11/// eigenvalue problems. The class provides methods to perform the factorization,
12/// check if the matrix is positive definite, calculate a witness vector if it is
13/// not positive definite, and calculate the symmetric quadratic form.
14///
15/// - LDL^T square-root-free version
16/// - Option allow semidefinite
17/// - A matrix A in R^{m x m} is positive definite iff v' A v > 0
18/// for all v in R^n.
19/// - O(p^2) per iteration, independent of N
20///
21/// Properties:
22///
23/// * `pos`: A tuple containing two usize values. It represents the dimensions of the LDLTMgr factorization.
24/// * `wit`: The `wit` property is an Array1 of f64 values.
25/// * `ndim`: The `ndim` property represents the size of the matrix that will be factorized using LDLTMgr
26/// factorization.
27/// * `storage`: The `storage` property is a 2-dimensional array of type `f64`. It is used to store the LDLTMgr
28/// factorization of a matrix.
29pub struct LDLTMgr {
30 pub pos: (usize, usize),
31 pub wit: Array1<f64>,
32 pub ndim: usize,
33 pub storage: Array2<f64>,
34}
35
36impl LDLTMgr {
37 /// The function `new` initializes a struct with default values for its fields.
38 ///
39 /// Arguments:
40 ///
41 /// * `ndim`: The parameter `ndim` represents the size of the arrays and matrices being initialized in the
42 /// `new` function. It is of type `usize`, which means it represents a non-negative integer.
43 ///
44 /// Returns:
45 ///
46 /// The `new` function is returning an instance of the struct that it is defined in.
47 ///
48 /// # Examples
49 ///
50 /// ```
51 /// use ellalgo_rs::oracles::ldlt_mgr::LDLTMgr;
52 /// use ndarray::{array, Array1};
53 /// let ldlt_mgr = LDLTMgr::new(3);
54 /// assert_eq!(ldlt_mgr.pos, (0, 0));
55 /// assert_eq!(ldlt_mgr.wit.len(), 3);
56 /// assert_eq!(ldlt_mgr.ndim, 3);
57 /// assert_eq!(ldlt_mgr.storage.len(), 9);
58 /// ```
59 pub fn new(ndim: usize) -> Self {
60 Self {
61 pos: (0, 0),
62 wit: Array1::zeros(ndim),
63 ndim,
64 storage: Array2::zeros((ndim, ndim)),
65 }
66 }
67
68 /// The `factorize` function takes a 2D array of f64 values and factors it using a closure.
69 ///
70 /// Arguments:
71 ///
72 /// * `mat_a`: The parameter `mat_a` is a reference to a 2-dimensional array (`Array2`) of `f64` values.
73 ///
74 /// Returns:
75 ///
76 /// The `factorize` function returns a boolean value.
77 ///
78 /// # Examples
79 ///
80 /// ```
81 /// use ellalgo_rs::oracles::ldlt_mgr::LDLTMgr;
82 /// use ndarray::array;
83 /// let mut ldlt_mgr = LDLTMgr::new(3);
84 /// let mat_a = array![
85 /// [1.0, 2.0, 3.0],
86 /// [2.0, 4.0, 5.0],
87 /// [3.0, 5.0, 6.0],
88 /// ];
89 /// assert_eq!(ldlt_mgr.factorize(&mat_a), false);
90 /// ```
91 pub fn factorize(&mut self, mat_a: &Array2<f64>) -> bool {
92 self.factor(|i, j| mat_a[[i, j]])
93 }
94
95 /// The `factor` function performs LDLTMgr factorization on a matrix and checks if it is symmetric
96 /// positive definite.
97 ///
98 /// Arguments:
99 ///
100 /// * `get_elem`: `get_elem` is a closure that takes two `usize` parameters (`i` and `j`) and returns a
101 /// `f64` value. It is used to retrieve elements from a matrix-like data structure.
102 ///
103 /// Returns:
104 ///
105 /// The `factor` function returns a boolean value.
106 ///
107 /// # Examples
108 ///
109 /// ```
110 /// use ellalgo_rs::oracles::ldlt_mgr::LDLTMgr;
111 /// use ndarray::array;
112 /// let mut ldlt_mgr = LDLTMgr::new(3);
113 /// let mat_a = array![
114 /// [1.0, 2.0, 3.0],
115 /// [2.0, 4.0, 5.0],
116 /// [3.0, 5.0, 6.0],
117 /// ];
118 /// assert_eq!(ldlt_mgr.factor(&|i, j| mat_a[[i, j]]), false);
119 /// ```
120 pub fn factor<F>(&mut self, get_elem: F) -> bool
121 where
122 F: Fn(usize, usize) -> f64,
123 {
124 self.pos = (0, 0);
125 for i in 0..self.ndim {
126 let mut diag = get_elem(i, 0);
127 for j in 0..i {
128 self.storage[[j, i]] = diag;
129 self.storage[[i, j]] = diag / self.storage[[j, j]];
130 let stop = j + 1;
131 // diag = get_elem(i, stop);
132 // for k in 0..stop {
133 // diag -= self.storage[[i, k]] * self.storage[[k, stop]];
134 // }
135 diag = get_elem(i, stop)
136 - self
137 .storage
138 .slice(s![i, 0..stop])
139 .dot(&self.storage.slice(s![0..stop, stop]));
140 }
141 self.storage[[i, i]] = diag;
142 if diag <= 0.0 {
143 self.pos = (0, i + 1);
144 break;
145 }
146 }
147 self.is_spd()
148 }
149
150 /// The function `factor_with_allow_semidefinite` checks if a given matrix is symmetric positive
151 /// definite (SPD) or semidefinite.
152 ///
153 /// Arguments:
154 ///
155 /// * `get_elem`: `get_elem` is a closure that takes two `usize` parameters `i` and `start` and returns
156 /// a `f64` value.
157 ///
158 /// Returns:
159 ///
160 /// a boolean value.
161 ///
162 /// # Examples
163 ///
164 /// ```
165 /// use ellalgo_rs::oracles::ldlt_mgr::LDLTMgr;
166 /// use ndarray::array;
167 /// let mut ldlt_mgr = LDLTMgr::new(3);
168 /// let mat_a = array![
169 /// [1.0, 2.0, 3.0],
170 /// [2.0, 4.0, 5.0],
171 /// [3.0, 5.0, 6.0],
172 /// ];
173 /// assert_eq!(ldlt_mgr.factor_with_allow_semidefinite(&|i, j| mat_a[[i, j]]), true);
174 /// ```
175 pub fn factor_with_allow_semidefinite<F>(&mut self, get_elem: F) -> bool
176 where
177 F: Fn(usize, usize) -> f64,
178 {
179 self.pos = (0, 0);
180 let mut start = 0;
181 for i in 0..self.ndim {
182 let mut diag = get_elem(i, start);
183 for j in start..i {
184 self.storage[[j, i]] = diag;
185 self.storage[[i, j]] = diag / self.storage[[j, j]];
186 let stop = j + 1;
187 // diag = get_elem(i, stop);
188 // for k in start..stop {
189 // diag -= self.storage[[i, k]] * self.storage[[k, stop]];
190 // }
191 diag = get_elem(i, stop)
192 - self
193 .storage
194 .slice(s![i, start..stop])
195 .dot(&self.storage.slice(s![start..stop, stop]));
196 }
197 self.storage[[i, i]] = diag;
198 if diag < 0.0 {
199 self.pos = (start, i + 1);
200 break;
201 }
202 if diag == 0.0 {
203 start = i + 1;
204 // restart at i + 1, special as an LMI oracle
205 }
206 }
207 self.is_spd()
208 }
209
210 /// The function `is_spd` checks if the second element of the tuple `pos` is equal to 0.
211 ///
212 /// Returns:
213 ///
214 /// A boolean value is being returned.
215 ///
216 /// # Examples
217 ///
218 /// ```
219 /// use ellalgo_rs::oracles::ldlt_mgr::LDLTMgr;
220 /// use ndarray::array;
221 /// let mut ldlt_mgr = LDLTMgr::new(3);
222 /// let mat_a = array![
223 /// [1.0, 2.0, 3.0],
224 /// [2.0, 4.0, 5.0],
225 /// [3.0, 5.0, 6.0],
226 /// ];
227 /// assert_eq!(ldlt_mgr.factor(&|i, j| mat_a[[i, j]]), false);
228 /// assert_eq!(ldlt_mgr.is_spd(), false);
229 /// ```
230 pub fn is_spd(&self) -> bool {
231 self.pos.1 == 0
232 }
233
234 /// The function calculates the witness vector of a matrix.
235 ///
236 /// Returns:
237 ///
238 /// The function `witness` returns a `f64` (a 64-bit floating-point number).
239 ///
240 /// # Examples
241 ///
242 /// ```
243 /// use ellalgo_rs::oracles::ldlt_mgr::LDLTMgr;
244 /// use ndarray::array;
245 /// let mut ldlt_mgr = LDLTMgr::new(3);
246 /// let mat_a = array![
247 /// [1.0, 2.0, 3.0],
248 /// [2.0, 4.0, 5.0],
249 /// [3.0, 5.0, 6.0],
250 /// ];
251 /// assert_eq!(ldlt_mgr.factor(&|i, j| mat_a[[i, j]]), false);
252 /// assert_eq!(ldlt_mgr.witness(), 0.0);
253 /// assert_eq!(ldlt_mgr.wit[0], -2.0);
254 /// assert_eq!(ldlt_mgr.wit[1], 1.0);
255 /// assert_eq!(ldlt_mgr.pos.1, 2);
256 /// ```
257 pub fn witness(&mut self) -> f64 {
258 if self.is_spd() {
259 panic!("Matrix is symmetric positive definite");
260 }
261 let (start, ndim) = self.pos;
262 let m = ndim - 1;
263 self.wit[m] = 1.0;
264 for i in (start..m).rev() {
265 self.wit[i] = 0.0;
266 for k in i..ndim {
267 self.wit[i] -= self.storage[[k, i]] * self.wit[k];
268 }
269 }
270 -self.storage[[m, m]]
271 }
272
273 /// The `sym_quad` function calculates the quadratic form of a symmetric matrix and a vector.
274 ///
275 /// Arguments:
276 ///
277 /// * `mat_a`: A 2-dimensional array of type f64.
278 ///
279 /// Returns:
280 ///
281 /// The function `sym_quad` returns a `f64` value.
282 ///
283 /// # Examples
284 ///
285 /// ```
286 /// use ellalgo_rs::oracles::ldlt_mgr::LDLTMgr;
287 /// use ndarray::array;
288 /// let mut ldlt_mgr = LDLTMgr::new(3);
289 /// let mat_a = array![
290 /// [1.0, 2.0, 3.0],
291 /// [2.0, 4.0, 5.0],
292 /// [3.0, 5.0, 6.0],
293 /// ];
294 /// assert_eq!(ldlt_mgr.factor(&|i, j| mat_a[[i, j]]), false);
295 /// assert_eq!(ldlt_mgr.witness(), 0.0);
296 /// let mat_b = array![
297 /// [1.0, 2.0, 3.0],
298 /// [2.0, 6.0, 5.0],
299 /// [3.0, 5.0, 4.0],
300 /// ];
301 /// assert_eq!(ldlt_mgr.sym_quad(&mat_b), 2.0);
302 pub fn sym_quad(&self, mat_a: &Array2<f64>) -> f64 {
303 let mut res = 0.0;
304 let (start, stop) = self.pos;
305 for i in start..stop {
306 let mut s = 0.0;
307 for j in (i + 1)..stop {
308 s += mat_a[[i, j]] * self.wit[j];
309 }
310 res += self.wit[i] * (mat_a[[i, i]] * self.wit[i] + 2.0 * s);
311 }
312 res
313 }
314
315 pub fn sqrt(&self) -> Array2<f64> {
316 if !self.is_spd() {
317 panic!("Matrix is not symmetric positive definite");
318 }
319 let mut res = Array2::zeros((self.ndim, self.ndim));
320 for i in 0..self.ndim {
321 res[[i, i]] = self.storage[[i, i]].sqrt();
322 for j in (i + 1)..self.ndim {
323 res[[i, j]] = self.storage[[j, i]] * res[[i, i]];
324 }
325 }
326 res
327 }
328}
329
330// pub fn test_ldlt() {
331// let ndim = 3;
332// let mut ldlt_mgr = LDLTMgr::new(ndim);
333// let a = Array::from_shape_vec((ndim, ndim), vec![4.0, 12.0, -16.0, 12.0, 37.0, -43.0, -16.0, -43.0, 98.0]).unwrap();
334// ldlt_mgr.factorize(&a);
335// assert!(ldlt_mgr.is_spd());
336// let b = Array::from_shape_vec((ndim, ndim), vec![4.0, 12.0, -16.0, 12.0, 37.0, -43.0, -16.0, -43.0, 99.0]).unwrap();
337// ldlt_mgr.factorize(&b);
338// assert!(!ldlt_mgr.is_spd());
339// assert_eq!(ldlt_mgr.witness(), -1.0);
340// assert_eq!(ldlt_mgr.sym_quad(&a), 16.0);
341// }
342
343#[cfg(test)]
344mod tests {
345 use super::*;
346 use ndarray::{Array2, ShapeError};
347
348 #[test]
349 fn test_chol1() -> Result<(), ShapeError> {
350 let l1 = Array2::from_shape_vec(
351 (3, 3),
352 vec![25.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0],
353 )?;
354 let mut ldlt_mgr = LDLTMgr::new(3);
355 assert!(ldlt_mgr.factorize(&l1));
356 Ok(())
357 }
358
359 #[test]
360 fn test_chol2() -> Result<(), ShapeError> {
361 let l2 = Array2::from_shape_vec(
362 (4, 4),
363 vec![
364 18.0, 22.0, 54.0, 42.0, 22.0, -70.0, 86.0, 62.0, 54.0, 86.0, -174.0, 134.0, 42.0,
365 62.0, 134.0, -106.0,
366 ],
367 )?;
368 let mut ldlt_mgr = LDLTMgr::new(4);
369 assert!(!ldlt_mgr.factorize(&l2));
370 ldlt_mgr.witness();
371 assert_eq!(ldlt_mgr.pos, (0, 2));
372 Ok(())
373 }
374
375 #[test]
376 fn test_chol3() -> Result<(), ShapeError> {
377 let l3 = Array2::from_shape_vec(
378 (3, 3),
379 vec![0.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0],
380 )?;
381 let mut ldlt_mgr = LDLTMgr::new(3);
382 assert!(!ldlt_mgr.factorize(&l3));
383 let ep = ldlt_mgr.witness();
384 assert_eq!(ldlt_mgr.pos, (0, 1));
385 assert_eq!(ldlt_mgr.wit[0], 1.0);
386 assert_eq!(ep, 0.0);
387 Ok(())
388 }
389
390 #[test]
391 fn test_chol4() -> Result<(), ShapeError> {
392 let l1 = Array2::from_shape_vec(
393 (3, 3),
394 vec![25.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0],
395 )?;
396 let mut ldlt_mgr = LDLTMgr::new(3);
397 assert!(ldlt_mgr.factor_with_allow_semidefinite(|i, j| l1[[i, j]]));
398 Ok(())
399 }
400
401 #[test]
402 fn test_chol5() -> Result<(), ShapeError> {
403 let l2 = Array2::from_shape_vec(
404 (4, 4),
405 vec![
406 18.0, 22.0, 54.0, 42.0, 22.0, -70.0, 86.0, 62.0, 54.0, 86.0, -174.0, 134.0, 42.0,
407 62.0, 134.0, -106.0,
408 ],
409 )?;
410 let mut ldlt_mgr = LDLTMgr::new(4);
411 assert!(!ldlt_mgr.factor_with_allow_semidefinite(|i, j| l2[[i, j]]));
412 ldlt_mgr.witness();
413 assert_eq!(ldlt_mgr.pos, (0, 2));
414 Ok(())
415 }
416
417 #[test]
418 fn test_chol6() -> Result<(), ShapeError> {
419 let l3 = Array2::from_shape_vec(
420 (3, 3),
421 vec![0.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0],
422 )?;
423 let mut ldlt_mgr = LDLTMgr::new(3);
424 assert!(ldlt_mgr.factor_with_allow_semidefinite(|i, j| l3[[i, j]]));
425 Ok(())
426 }
427
428 #[test]
429 fn test_chol7() -> Result<(), ShapeError> {
430 let l3 = Array2::from_shape_vec(
431 (3, 3),
432 vec![0.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, -20.0],
433 )?;
434 let mut ldlt_mgr = LDLTMgr::new(3);
435 assert!(!ldlt_mgr.factor_with_allow_semidefinite(|i, j| l3[[i, j]]));
436 let ep = ldlt_mgr.witness();
437 assert_eq!(ep, 20.0);
438 Ok(())
439 }
440
441 #[test]
442 fn test_chol8() -> Result<(), ShapeError> {
443 let l3 = Array2::from_shape_vec(
444 (3, 3),
445 vec![0.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 20.0],
446 )?;
447 let mut ldlt_mgr = LDLTMgr::new(3);
448 assert!(!ldlt_mgr.factorize(&l3));
449 Ok(())
450 }
451
452 #[test]
453 fn test_chol9() -> Result<(), ShapeError> {
454 let l3 = Array2::from_shape_vec(
455 (3, 3),
456 vec![0.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 20.0],
457 )?;
458 let mut ldlt_mgr = LDLTMgr::new(3);
459 assert!(ldlt_mgr.factor_with_allow_semidefinite(|i, j| l3[[i, j]]));
460 Ok(())
461 }
462
463 #[test]
464 fn test_ldlt_mgr_sqrt() -> Result<(), ShapeError> {
465 let a =
466 Array2::from_shape_vec((3, 3), vec![1.0, 0.5, 0.5, 0.5, 1.25, 0.75, 0.5, 0.75, 1.5])?;
467 let mut ldlt_mgr = LDLTMgr::new(3);
468 ldlt_mgr.factorize(&a);
469 assert!(ldlt_mgr.is_spd());
470 let r = ldlt_mgr.sqrt();
471 assert_eq!(
472 r,
473 Array2::from_shape_vec((3, 3), vec![1.0, 0.5, 0.5, 0.0, 1.0, 0.5, 0.0, 0.0, 1.0])?
474 );
475 Ok(())
476 }
477}