jf_pcs/
toeplitz.rs

1// Copyright (c) 2023 Espresso Systems (espressosys.com)
2// This file is part of the Jellyfish library.
3
4// You should have received a copy of the MIT License
5// along with the Jellyfish library. If not, see <https://mit-license.org/>.
6
7//! Toeplitz matrices and Circulant matrices operations
8//! References: `https://eprint.iacr.org/2020/1516.pdf`
9
10use ark_ff::FftField;
11use ark_poly::{domain::DomainCoeff, EvaluationDomain, GeneralEvaluationDomain};
12use ark_std::{format, ops::Mul, string::ToString, vec::Vec};
13use jf_utils::hadamard_product;
14
15use crate::prelude::PCSError;
16
17/// An `NxN` [Circulant Matrix](https://en.wikipedia.org/wiki/Circulant_matrix)
18/// is unambiguously represented by its first column, and has the form:
19///
20/// [c_0   c_N-1 c_N-2 ... c_1]
21/// [c_1   c_0   c_N-1 ... c_2]
22/// [c_2   c_1   c_0   ... c_3]
23/// [..    ..    ..    ... .. ]
24/// [c_N-1 c_N-2 c_N-3 ... c_0]
25///
26/// It is a special case of [`ToeplitzMatrix`].
27#[derive(Debug, Clone, PartialEq, Eq)]
28pub struct CirculantMatrix<F: FftField> {
29    col: Vec<F>,
30}
31
32impl<F: FftField> CirculantMatrix<F> {
33    /// Construct a Circulant matrix by its first column vector
34    pub fn new(col: Vec<F>) -> Self {
35        Self { col }
36    }
37
38    /// Fast multiplication of a Circulant matrix by a vector via FFT
39    /// Details see Section 2.2.1 of [Tomescu20](https://eprint.iacr.org/2020/1516.pdf).
40    // TODO: (alex) think of ways to extend to arbitrary length vector (simply
41    // truncate doesn't work).
42    pub fn fast_vec_mul<T>(&self, m: &[T]) -> Result<Vec<T>, PCSError>
43    where
44        T: for<'a> Mul<&'a F, Output = T> + DomainCoeff<F>,
45    {
46        if !m.len().is_power_of_two() {
47            return Err(PCSError::InvalidParameters(
48                "Fast Circulant Matrix mul only supports vector size of power of two.".to_string(),
49            ));
50        }
51        if m.len() != self.col.len() {
52            return Err(PCSError::InvalidParameters(
53                "Wrong input dimension for matrix mul.".to_string(),
54            ));
55        }
56        let domain: GeneralEvaluationDomain<F> =
57            GeneralEvaluationDomain::new(m.len()).expect("Should init an evaluation domain");
58
59        let m_evals = domain.fft(m); // DFT(m)
60        let col_evals = domain.fft(&self.col); // DFT(c_N)
61        let eval_prod = // DFT(c_N) * DFT(m)
62            hadamard_product(col_evals, m_evals).expect("Hadamard product should succeed");
63        let res = domain.ifft(&eval_prod); // iDFT(DFT(c_N) * DFT(m))
64        Ok(res)
65    }
66}
67
68/// An `NxN` [Toeplitz Matrix](https://en.wikipedia.org/wiki/Toeplitz_matrix) is
69/// unambiguously represented by its first column and first row, and has the
70/// form:
71///
72/// [a_0   a_-1  a_-2  ...  a_-(N-1)]
73/// [a_1   a_0   a_-1  ...  a_-(N-2)]
74/// [a_2   a_1   a_0   ...  a_-(N-3)]
75/// [..    ..    ..    ...    ..    ]
76/// [a_N-1 a_N-2 ..    ...  a_0     ]
77#[derive(Debug, Clone, PartialEq, Eq)]
78pub struct ToeplitzMatrix<F: FftField> {
79    col: Vec<F>,
80    row: Vec<F>,
81}
82
83impl<F: FftField> ToeplitzMatrix<F> {
84    /// constructor for a new Toeplitz matrix.
85    pub fn new(col: Vec<F>, row: Vec<F>) -> Result<Self, PCSError> {
86        if col.is_empty() || col.len() != row.len() {
87            return Err(PCSError::InvalidParameters(format!(
88                "row: {}, col: {} should be both positive and equal",
89                row.len(),
90                col.len()
91            )));
92        }
93        if col[0] != row[0] {
94            return Err(PCSError::InvalidParameters(format!(
95                "1st value in 1st col: {:?} should be the same as that in 1st row: {:?}",
96                col[0], row[0]
97            )));
98        }
99
100        Ok(Self { col, row })
101    }
102
103    /// Embeds a Toeplitz matrix of size N to a Circulant matrix of size 2N.
104    ///
105    /// Details see Section 2.3.1 of [Tomescu20](https://eprint.iacr.org/2020/1516.pdf).
106    // TODO: (alex) turn this into concurrent code after: https://github.com/EspressoSystems/jellyfish/issues/111
107    pub fn circulant_embedding(&self) -> Result<CirculantMatrix<F>, PCSError> {
108        let mut extension_col = self.row.clone();
109        extension_col.rotate_left(1);
110        extension_col.reverse();
111
112        Ok(CirculantMatrix {
113            col: [self.col.clone(), extension_col].concat(),
114        })
115    }
116
117    /// Fast multiplication of a Toeplitz matrix by embedding it into a
118    /// circulant matrix and multiply there.
119    ///
120    /// Details see Section 2.3.1 of [Tomescu20](https://eprint.iacr.org/2020/1516.pdf).
121    pub fn fast_vec_mul<T>(&self, v: &[T]) -> Result<Vec<T>, PCSError>
122    where
123        T: for<'a> Mul<&'a F, Output = T> + DomainCoeff<F>,
124    {
125        if !v.len().is_power_of_two() {
126            return Err(PCSError::InvalidParameters(
127                "Fast Toeplitz Matrix mul only supports vector size of power of two.".to_string(),
128            ));
129        }
130        if v.len() != self.col.len() {
131            return Err(PCSError::InvalidParameters(
132                "Wrong input dimension for matrix mul.".to_string(),
133            ));
134        }
135
136        let cir_repr = self.circulant_embedding()?;
137        let mut padded_v = Vec::from(v);
138        padded_v.resize(2 * v.len(), T::zero());
139
140        let mut res = cir_repr.fast_vec_mul(&padded_v)?;
141        res.truncate(v.len());
142        Ok(res)
143    }
144}
145
146impl<F: FftField> From<CirculantMatrix<F>> for ToeplitzMatrix<F> {
147    fn from(c: CirculantMatrix<F>) -> Self {
148        let mut row = c.col.clone();
149        row.rotate_left(1);
150        row.reverse();
151
152        Self { col: c.col, row }
153    }
154}
155
156impl<F: FftField> TryFrom<ToeplitzMatrix<F>> for CirculantMatrix<F> {
157    type Error = PCSError;
158
159    fn try_from(t: ToeplitzMatrix<F>) -> Result<Self, Self::Error> {
160        let mut expected_col = t.row;
161        expected_col.reverse();
162        expected_col.rotate_right(1);
163        if expected_col != t.col {
164            return Err(PCSError::InvalidParameters(
165                "Not a Circulant Matrix".to_string(),
166            ));
167        }
168        Ok(Self { col: t.col })
169    }
170}
171
172#[cfg(test)]
173mod tests {
174    use super::*;
175    use ark_bls12_381::{Fr, G1Projective};
176    use ark_ff::Field;
177    use ark_std::{ops::AddAssign, vec, UniformRand};
178    use jf_utils::test_rng;
179
180    // a MxN matrix, M rows, N cols.
181    struct Matrix<T, const M: usize, const N: usize>([[T; N]; M]);
182    struct Vector<T, const N: usize>([T; N]);
183
184    impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
185        fn transpose(self) -> Matrix<T, N, M> {
186            let mut transposed = [[self.0[0][0]; M]; N];
187
188            #[allow(clippy::needless_range_loop)]
189            for i in 0..M {
190                for j in 0..N {
191                    transposed[j][i] = self.0[i][j];
192                }
193            }
194            Matrix(transposed)
195        }
196    }
197
198    impl<T: Copy, const N: usize> From<Matrix<T, 1, N>> for Vector<T, N> {
199        fn from(m: Matrix<T, 1, N>) -> Self {
200            let mut v = [m.0[0][0]; N];
201            v.copy_from_slice(&m.0[0]);
202            Vector(v)
203        }
204    }
205
206    impl<T: Copy, const N: usize> From<Vector<T, N>> for Vec<T> {
207        fn from(v: Vector<T, N>) -> Self {
208            Vec::from(&v.0[..])
209        }
210    }
211
212    impl<F: FftField> CirculantMatrix<F> {
213        fn full_matrix<const N: usize>(self) -> Matrix<F, N, N> {
214            assert_eq!(self.col.len(), N);
215            let t: ToeplitzMatrix<F> = self.into();
216
217            let mut first_row = [F::default(); N];
218            first_row.copy_from_slice(&t.row);
219            let mut row_vecs = [first_row; N];
220            let mut cur_row = first_row;
221
222            #[allow(clippy::needless_range_loop)]
223            for i in 1..N {
224                cur_row.rotate_right(1);
225                row_vecs[i] = cur_row;
226            }
227            // some arbitrary sanity check
228            assert_eq!(row_vecs[N - 1][0], row_vecs[0][1]);
229            assert_eq!(row_vecs[1][0], row_vecs[0][N - 1]);
230            assert_eq!(row_vecs[N - 1][N - 1], row_vecs[0][0]);
231
232            Matrix(row_vecs)
233        }
234    }
235
236    impl<F: FftField> ToeplitzMatrix<F> {
237        fn full_matrix<const N: usize>(self) -> Matrix<F, N, N> {
238            assert_eq!(self.col.len(), N);
239            let mut matrix = [[F::zero(); N]; N];
240            for i in 0..N {
241                matrix[i][0] = self.col[i];
242                matrix[0][i] = self.row[i];
243            }
244            for i in 1..N {
245                for j in 1..N {
246                    matrix[i][j] = matrix[i - 1][j - 1];
247                }
248            }
249            Matrix(matrix)
250        }
251    }
252
253    fn naive_matrix_mul<F, T, const M: usize, const N: usize, const K: usize>(
254        a: Matrix<F, M, N>,
255        b: Matrix<T, N, K>,
256    ) -> Matrix<T, M, K>
257    where
258        F: Field,
259        T: for<'a> Mul<&'a F, Output = T> + AddAssign<T> + Copy + Default,
260    {
261        let mut c = [[T::default(); K]; M];
262
263        #[allow(clippy::needless_range_loop)]
264        for i in 0..M {
265            for j in 0..K {
266                for k in 0..N {
267                    c[i][j] += b.0[k][j] * &a.0[i][k];
268                }
269            }
270        }
271        Matrix(c)
272    }
273
274    #[test]
275    fn test_circulant_mul() -> Result<(), PCSError> {
276        let mut rng = test_rng();
277        // happy path
278        const N: usize = 16;
279
280        let cir_matrix = CirculantMatrix::new((0..N).map(|_| Fr::rand(&mut rng)).collect());
281
282        let msgs = [G1Projective::rand(&mut rng); N];
283        let msg_matrix = Matrix([msgs]);
284
285        let expected: Vector<G1Projective, N> =
286            naive_matrix_mul(cir_matrix.clone().full_matrix(), msg_matrix.transpose())
287                .transpose()
288                .into();
289        let got = cir_matrix.fast_vec_mul(&msgs)?;
290        assert_eq!(
291            <Vector<G1Projective, N> as Into<Vec<G1Projective>>>::into(expected),
292            got,
293            "Fast Circulant Matrix mul for EC group is incorrect."
294        );
295
296        let f_msgs = [Fr::rand(&mut rng); N];
297        let f_msg_matrix = Matrix([f_msgs]);
298
299        let expected: Vector<Fr, N> =
300            naive_matrix_mul(cir_matrix.clone().full_matrix(), f_msg_matrix.transpose())
301                .transpose()
302                .into();
303        let got = cir_matrix.fast_vec_mul(&f_msgs)?;
304        assert_eq!(
305            <Vector<Fr, N> as Into<Vec<Fr>>>::into(expected),
306            got,
307            "Fast Circulant Matrix mul for field is incorrect."
308        );
309
310        // bad path
311        // mismatched matrix.col.len() and msgs.len() should fail
312        let bad_msg = [msgs.to_vec(), vec![G1Projective::rand(&mut rng)]].concat();
313        assert!(cir_matrix.fast_vec_mul(&bad_msg).is_err());
314
315        // non power-of-two matrix fast mul should fail
316        let m = bad_msg.len(); // same dimension as the message, but not a power-of-two
317        let cir_matrix = CirculantMatrix::new((0..m).map(|_| Fr::rand(&mut rng)).collect());
318
319        assert!(
320            !m.is_power_of_two()
321                && m == cir_matrix.col.len()
322                && cir_matrix.fast_vec_mul(&bad_msg).is_err()
323        );
324
325        Ok(())
326    }
327
328    #[test]
329    fn test_toeplitz_mul() -> Result<(), PCSError> {
330        let mut rng = test_rng();
331        const N: usize = 16;
332
333        let rand_col: Vec<Fr> = (0..N).map(|_| Fr::rand(&mut rng)).collect();
334        let rand_row = (0..N)
335            .map(|i| {
336                if i == 0 {
337                    rand_col[0]
338                } else {
339                    Fr::rand(&mut rng)
340                }
341            })
342            .collect();
343        let toep_matrix = ToeplitzMatrix::new(rand_col, rand_row)?;
344
345        let mut msgs = [G1Projective::default(); N];
346        for m in msgs.iter_mut() {
347            *m = G1Projective::rand(&mut rng);
348        }
349        let msg_matrix = Matrix([msgs]);
350
351        let expected: Vector<G1Projective, N> =
352            naive_matrix_mul(toep_matrix.clone().full_matrix(), msg_matrix.transpose())
353                .transpose()
354                .into();
355        let got = toep_matrix.fast_vec_mul(&msgs)?;
356        assert_eq!(
357            <Vector<G1Projective, N> as Into<Vec<G1Projective>>>::into(expected),
358            got,
359            "Fast Toeplitz Matrix mul is incorrect."
360        );
361
362        Ok(())
363    }
364}