optimization_engine/
matrix_operations.rs

1//! # matrix_operations
2//!
3//! Simple algebraic operations used by Optimization Engine
4//!
5//! # Examples
6//!
7//! ```
8//! use optimization_engine::matrix_operations::*;
9//!
10//! let a = [1.0, 2.0, 3.0];
11//! let b = [4.0, 5.0, 6.0];
12//!
13//! // inner product
14//! let a_dot_b = inner_product(&a, &b);
15//! assert!(a_dot_b == 32.);
16//!
17//! // Euclidean norm and squared norm
18//! let norm_a = norm2(&a);
19//! let norm_sq_a = norm2_squared(&a);
20//! assert!(norm_sq_a == 14.);
21//!
22//! // Squared Euclidean norm of difference of vectors
23//! let norm_sq_a_minus_b = norm2_squared_diff(&a, &b);
24//! assert!(norm_sq_a_minus_b == 27.);
25//!
26//! // Sum of elements of vector
27//! let sum_a = sum(&a);
28//! assert!(sum_a == 6.);
29//!
30//! // Check whether all elements of a vector are finite
31//! assert!(is_finite(&a));
32//!
33//! // Infinity norm
34//! let norm_inf_b = norm_inf(&b);
35//! assert!(norm_inf_b == 6.);
36//! ```
37//!
38
39use num::{Float, Zero};
40use std::iter::Sum;
41use std::ops::Mul;
42
43/// Calculate the inner product of two vectors
44#[inline(always)]
45pub fn inner_product<T>(a: &[T], b: &[T]) -> T
46where
47    T: Float + Sum<T> + Mul<T, Output = T>,
48{
49    assert!(a.len() == b.len());
50
51    a.iter().zip(b.iter()).map(|(x, y)| (*x) * (*y)).sum()
52}
53
54/// Calculate the 1-norm of a vector
55#[inline(always)]
56pub fn norm1<T>(a: &[T]) -> T
57where
58    T: Float + Sum<T>,
59{
60    a.iter().map(|x| x.abs()).sum()
61}
62
63/// Calculate the 2-norm of a vector
64#[inline(always)]
65pub fn norm2<T>(a: &[T]) -> T
66where
67    T: Float + Sum<T> + Mul<T, Output = T>,
68{
69    let norm: T = norm2_squared(a);
70    norm.sqrt()
71}
72
73/// Calculate the squared 2-norm of the difference of two vectors
74#[inline(always)]
75pub fn norm2_squared_diff<T>(a: &[T], b: &[T]) -> T
76where
77    T: Float + Sum<T> + Mul<T, Output = T> + std::ops::AddAssign,
78{
79    a.iter().zip(b.iter()).fold(T::zero(), |mut sum, (&x, &y)| {
80        sum += (x - y).powi(2);
81        sum
82    })
83}
84
85/// Calculate the 2-norm of a vector
86#[inline(always)]
87pub fn norm2_squared<T>(a: &[T]) -> T
88where
89    T: Float + Sum<T> + Mul<T, Output = T>,
90{
91    let norm: T = a.iter().map(|x| (*x) * (*x)).sum();
92    norm
93}
94
95/// Calculate the sum of all elements of a vector
96#[inline(always)]
97pub fn sum<T>(a: &[T]) -> T
98where
99    T: Float + Sum<T> + Mul<T, Output = T>,
100{
101    let norm: T = a.iter().copied().sum();
102    norm
103}
104
105/// Calculates the infinity-norm of a vector
106#[inline(always)]
107pub fn norm_inf<T>(a: &[T]) -> T
108where
109    T: Float + Zero,
110{
111    a.iter()
112        .fold(T::zero(), |current_max, x| x.abs().max(current_max))
113}
114
115/// Computes the infinity norm of the difference of two vectors
116#[inline(always)]
117pub fn norm_inf_diff<T>(a: &[T], b: &[T]) -> T
118where
119    T: Float + Zero,
120{
121    assert_eq!(a.len(), b.len());
122    a.iter()
123        .zip(b.iter())
124        .fold(T::zero(), |current_max, (x, y)| {
125            (*x - *y).abs().max(current_max)
126        })
127}
128
129/// Checks whether all elements of a vector are finite
130///
131/// ## Returns
132///
133/// Returns `true` if all elements are finite and `false` if any
134/// of the elements are either NaN or Infinity
135#[inline(always)]
136pub fn is_finite<T>(a: &[T]) -> bool
137where
138    T: Float,
139{
140    !a.iter().any(|&xi| !xi.is_finite())
141}
142
143/* ---------------------------------------------------------------------------- */
144/*          TESTS                                                               */
145/* ---------------------------------------------------------------------------- */
146#[cfg(test)]
147mod tests {
148    use crate::*;
149
150    #[test]
151    fn t_inner_product_test() {
152        unit_test_utils::assert_nearly_equal(
153            14.0,
154            matrix_operations::inner_product(&[1.0, 2.0, 3.0], &[1.0, 2.0, 3.0]),
155            1e-10,
156            1e-16,
157            "inner product",
158        );
159    }
160
161    #[test]
162    #[should_panic]
163    fn t_inner_product_test_panic() {
164        matrix_operations::inner_product(&[2.0, 3.0], &[1.0, 2.0, 3.0]);
165    }
166
167    #[test]
168    fn t_norm1_test() {
169        unit_test_utils::assert_nearly_equal(
170            6.0,
171            matrix_operations::norm1(&[1.0, -2.0, -3.0]),
172            1e-10,
173            1e-16,
174            "norm1",
175        );
176    }
177
178    #[test]
179    fn t_norm2_test() {
180        unit_test_utils::assert_nearly_equal(
181            5.0,
182            matrix_operations::norm2(&[3.0, 4.0]),
183            1e-10,
184            1e-16,
185            "norm2",
186        );
187    }
188
189    #[test]
190    fn t_norm_inf_test() {
191        unit_test_utils::assert_nearly_equal(
192            3.0,
193            matrix_operations::norm_inf(&[1.0, -2.0, -3.0]),
194            1e-10,
195            1e-16,
196            "norm infinity of vector",
197        );
198        unit_test_utils::assert_nearly_equal(
199            8.0,
200            matrix_operations::norm_inf(&[1.0, -8.0, -3.0, 0.0]),
201            1e-10,
202            1e-16,
203            "infinity norm",
204        );
205    }
206
207    #[test]
208    fn t_norm_inf_diff() {
209        let x = [1.0, 2.0, 1.0];
210        let y = [-4.0, 0.0, 3.0];
211        let norm_diff = matrix_operations::norm_inf_diff(&x, &y);
212        unit_test_utils::assert_nearly_equal(5.0f64, norm_diff, 1e-10, 1e-9, "norm of difference");
213        unit_test_utils::assert_nearly_equal(
214            0.0f64,
215            matrix_operations::norm_inf_diff(&x, &x),
216            1e-10,
217            1e-16,
218            "difference of same vector",
219        );
220        unit_test_utils::assert_nearly_equal(
221            0.0f64,
222            matrix_operations::norm_inf_diff(&[], &[]),
223            1e-10,
224            1e-16,
225            "difference of empty vectors",
226        );
227    }
228
229    #[test]
230    #[should_panic]
231    fn t_norm_inf_diff_panic() {
232        let x = [1.0, 2.0, 3.0];
233        let y = [0.0, 3.0];
234        let _ = matrix_operations::norm_inf_diff(&x, &y);
235    }
236
237    #[test]
238    fn t_norm2_squared_diff_test() {
239        let x = [2.0, 5.0, 7.0, -1.0];
240        let y = [4.0, 1.0, 0.0, 10.0];
241        let norm2sq = matrix_operations::norm2_squared_diff(&x, &y);
242        unit_test_utils::assert_nearly_equal(190., norm2sq, 1e-10, 1e-12, "norm sq diff");
243    }
244}