Marlin refactoring

README.md

@@ -31,11 +31,6 @@ The high-level structure of the repository is as follows:
   - `indexer.rs` , which contains structs and functions for the circuit-specific pre-processing stage. 
   - `prover.rs`  and `verifier.rs` , which implement the round functions of prover and verifier of the oracle proof.
 
-## Release Note ##
-
-The current release implements Coboundary Marlin using Marlin's original matrix arithmetization. 
-It is a refactoring of release 0.2.0.
-
 ## License
 
 This library is licensed under either of the following licenses, at your discretion.

src/iop/indexer.rs

@@ -97,15 +97,15 @@ pub struct Index<F: PrimeField> {
     /// The `C` matrix for the R1CS instance, in sparse representation
     pub c: SparseMatrix<F>,
 
-    /// Arithmetization of the kernel-matrix product `R_A`, which essentially contains
+    /// Arithmetization of the matrix`A`, which essentially contains
     /// the indexer polynomials `row(X)`, `col(X)`, `row.col(X)`, and `val.row.col(X)` of it.
-    pub a_star_arith: MatrixArithmetization<F>,
-    /// Arithmetization of the kernel-matrix product `R_B`, which essentially contains
+    pub a_arith: MatrixArithmetization<F>,
+    /// Arithmetization of the matrix`B`, which essentially contains
     /// the indexer polynomials `row(X)`, `col(X)`, `row.col(X)`, and `val.row.col(X)` of it.
-    pub b_star_arith: MatrixArithmetization<F>,
-    /// Arithmetization of the kernel-matrix product `R_C`, which essentially contains
+    pub b_arith: MatrixArithmetization<F>,
+    /// Arithmetization of the matrix`C`, which essentially contains
     /// the indexer polynomials `row(X)`, `col(X)`, `row.col(X)`, and `val.row.col(X)` of it.
-    pub c_star_arith: MatrixArithmetization<F>,
+    pub c_arith: MatrixArithmetization<F>,
 }
 
 impl<F: PrimeField> SemanticallyValid for Index<F> {
@@ -144,9 +144,7 @@ impl<F: PrimeField> SemanticallyValid for Index<F> {
                 && &m_star_arith.evals_on_domain_b.val_row_col.domain == &domain_b
         };
 
-        check_matrix(&self.a_star_arith)
-            && check_matrix(&self.b_star_arith)
-            && check_matrix(&self.c_star_arith)
+        check_matrix(&self.a_arith) && check_matrix(&self.b_arith) && check_matrix(&self.c_arith)
     }
 }
 
@@ -159,24 +157,28 @@ impl<F: PrimeField> Index<F> {
     /// Iterate over the indexed polynomials.
     pub fn iter(&self) -> impl Iterator<Item = &LabeledPolynomial<F>> {
         vec![
-            &self.a_star_arith.row,
-            &self.a_star_arith.col,
-            &self.a_star_arith.row_col,
-            &self.a_star_arith.val_row_col,
-            &self.b_star_arith.row,
-            &self.b_star_arith.col,
-            &self.b_star_arith.row_col,
-            &self.b_star_arith.val_row_col,
-            &self.c_star_arith.row,
-            &self.c_star_arith.col,
-            &self.c_star_arith.row_col,
-            &self.c_star_arith.val_row_col,
+            &self.a_arith.row,
+            &self.a_arith.col,
+            &self.a_arith.row_col,
+            &self.a_arith.val_row_col,
+            &self.b_arith.row,
+            &self.b_arith.col,
+            &self.b_arith.row_col,
+            &self.b_arith.val_row_col,
+            &self.c_arith.row,
+            &self.c_arith.col,
+            &self.c_arith.row_col,
+            &self.c_arith.val_row_col,
         ]
         .into_iter()
     }
 }
 
 impl<F: PrimeField> IOP<F> {
+    /// Build the four domains used in the protocol.
+    /// `num_aux` is the number of private witnesses
+    /// `num_inputs` is the number of public inputs (including the one for the constants)
+    /// `num_non_zero` is the max number of non-zero values in any of the matrices A, B, and C
     pub(crate) fn build_domains(
         num_aux: usize,
         num_inputs: usize,
@@ -217,7 +219,7 @@ impl<F: PrimeField> IOP<F> {
 
         // matrix post-processing: balance matrices
         let matrix_processing_time = start_timer!(|| "Processing matrices");
-        let (mut a, mut b, mut c) = post_process_matrices(&mut ics).expect("should not be `None`");
+        balance_matrices(&mut ics.at, &mut ics.bt);
         add_to_trace!(|| "number of (formatted) input_variables", || format!(
             "{}",
             ics.num_inputs
@@ -251,19 +253,19 @@ impl<F: PrimeField> IOP<F> {
             index_info.num_non_zero,
         )?;
 
+        let (mut a, mut b, mut c) =
+            post_process_matrices(&mut ics, &domain_h, &domain_x).expect("should not be `None`");
+
         let a_arithmetization_time = start_timer!(|| "Arithmetizing A");
-        let a_star_arith =
-            arithmetize_matrix("a", &mut a, &domain_k, &domain_h, &domain_x, &domain_b)?;
+        let a_star_arith = arithmetize_matrix("a", &mut a, &domain_k, &domain_h, &domain_b)?;
         end_timer!(a_arithmetization_time);
 
         let b_arithmetization_time = start_timer!(|| "Arithmetizing B");
-        let b_star_arith =
-            arithmetize_matrix("b", &mut b, &domain_k, &domain_h, &domain_x, &domain_b)?;
+        let b_star_arith = arithmetize_matrix("b", &mut b, &domain_k, &domain_h, &domain_b)?;
         end_timer!(b_arithmetization_time);
 
         let c_arithmetization_time = start_timer!(|| "Arithmetizing C");
-        let c_star_arith =
-            arithmetize_matrix("c", &mut c, &domain_k, &domain_h, &domain_x, &domain_b)?;
+        let c_star_arith = arithmetize_matrix("c", &mut c, &domain_k, &domain_h, &domain_b)?;
         end_timer!(c_arithmetization_time);
 
         end_timer!(index_time);
@@ -272,9 +274,9 @@ impl<F: PrimeField> IOP<F> {
             a,
             b,
             c,
-            a_star_arith,
-            b_star_arith,
-            c_star_arith,
+            a_arith: a_star_arith,
+            b_arith: b_star_arith,
+            c_arith: c_star_arith,
         })
     }
 }
@@ -285,11 +287,10 @@ impl<F: PrimeField> IOP<F> {
 
 ****************************************************************************/
 
-/// Contains information about the arithmetization of the kernel-matrix product
-/// `R_M(X,Y) = sum_{z in H} R(X,z) * M(z,Y)` of a sparse matrix `M`, as obtained
-/// by the lincheck to sumcheck reduction.
-/// The arithmetization is again with respect to the kernel `R(X,Y)`, i.e.
-///     `R_M(X,Y) = Sum_{z in K} val(z)*R(X,row(z))*R(Y,col(z))`
+/// Contains information about the arithmetization of a sparse matrix `M`,
+/// as obtained by the lincheck to sumcheck reduction.
+/// The arithmetization is with respect to the Lagrange kernel `L_H(X,Y)`, i.e.
+///     `M(X,Y) = Sum_{z in K} val(z) * L_H(X,row(z)) * L_H(Y,col(z))`
 /// over an *indexer domain* `K`, large enough to index the non-zero entries in
 /// `M`.
 #[derive(Derivative)]
@@ -345,16 +346,12 @@ pub struct MatrixEvals<F: PrimeField> {
     pub val_row_col: EvaluationsOnDomain<F>,
 }
 
-/// Given a sparse matrix `M`, computes the polynomial representations `val(X)`,
-/// `row(X)`, `col(X)`, and `row.col(X)` modulo `K` of the kernel representation of the
-/// derivative `R_M`. (Actually it returns the transpose, i.e. `row(X)` and `col(X)`
-/// interchanged)
-// The relation between the kernel representation and the Lagrange representation is in
-// general given by
-//     `M(X,Y) = sum_{x,y in H} m(x,y) * R(X,y)*R(x,Y) =
-//            = sum_{x,y in H} m(x,y)*R(x,x)*R(y,y) * L(X,y)*L(x,Y)`,
-// and hence `m(x,y)= M(x,y)/(R(x,x)R(y,y))` over `H x H`.
-// As `R_M(X,Y) = R(X,X)*M(X,Y)`, we overall have `m(x,y) = M(x,y)/R(y,y)`.
+/// Given a sparse matrix `M`, computes the polynomial representations `row(X)`, `col(X)`,
+/// `row.col(X)`, and `val.row.col(X)` of `M` such that
+///   M(X,Y) = sum_{w in K} val(w) * L_H(X,row(w)) * L_H(Y,col(w))
+/// where `K` is a domain large enough to index the non-zero entries of the matrix.
+/// In order to ease prover computations we provide `val.row.col(X)` instead of `val(X)`.
+/// For the same reason we also provide the polynomial `row.col(X)`.
 pub(crate) fn arithmetize_matrix<F: PrimeField>(
     matrix_name: &str,
     // The R1CS matrix.
@@ -363,54 +360,36 @@ pub(crate) fn arithmetize_matrix<F: PrimeField>(
     domain_k: &Box<dyn EvaluationDomain<F>>,
     // The domain `H` for the Lagrange representation of `M` .
     domain_h: &Box<dyn EvaluationDomain<F>>,
-    // The input domain `X`, a subdomain of the Lagrange domain `H`.
-    domain_x: &Box<dyn EvaluationDomain<F>>,
     // An extension of the indexer domain `K`, at least 4 times larger.
     domain_b: &Box<dyn EvaluationDomain<F>>,
 ) -> Result<MatrixArithmetization<F>, Error> {
-    let matrix_time = start_timer!(|| "Computing row, col, and val LDEs");
+    let matrix_time = start_timer!(|| format!("Arithemtizing matrix {}", matrix_name));
 
     let elems: Vec<_> = domain_h.elements().collect();
 
-    let mut row_vec = Vec::new();
-    let mut col_vec = Vec::new();
-    let mut val_vec = Vec::new();
+    let mut row_vec = Vec::with_capacity(domain_k.size());
+    let mut col_vec = Vec::with_capacity(domain_k.size());
+    let mut val_vec = Vec::with_capacity(domain_k.size());
 
-    let lde_evals_time = start_timer!(|| "Computing row, col and val evals");
+    let lde_evals_time = start_timer!(|| "Computing row, col, row.col and val.row.col evals");
 
-    let mut count = 0;
-
-    // As `R_M(X,Y) = R(X,X)*M(X,Y)`, we overall have `m(x,y) = M(x,y)/R(y,y)`.
     for (r, row) in matrix.into_iter().enumerate() {
         if !is_in_ascending_order(&row, |(_, a), (_, b)| a < b) {
             row.sort_by(|(_, a), (_, b)| a.cmp(b));
         };
 
-        for &mut (val, i) in row {
-            // As we do not re-index the y_A and y_B vectors by the input domain,
-            // we simply take elems[r]
-            let row_val = elems[r];
-            // on the contrary, column vectors are re-indexed
-            let col_val = elems[domain_h
-                .reindex_by_subdomain(domain_x.size(), i)
-                .map_err(|e| Error::Other(e.to_string()))?];
-
-            row_vec.push(row_val);
-            col_vec.push(col_val);
+        for &mut (val, c) in row {
+            row_vec.push(elems[r]);
+            col_vec.push(elems[c]);
             val_vec.push(val);
-
-            count += 1;
         }
     }
 
-    end_timer!(lde_evals_time);
+    // pad to len equal to domain_k.size()
+    row_vec.resize(domain_k.size(), elems[0]);
+    col_vec.resize(domain_k.size(), elems[0]);
+    val_vec.resize(domain_k.size(), F::zero());
 
-    // pad with zeroes
-    for _ in 0..(domain_k.size() - count) {
-        col_vec.push(elems[0]);
-        row_vec.push(elems[0]);
-        val_vec.push(F::zero());
-    }
     let (row_col_vec, val_row_col_vec): (Vec<_>, Vec<_>) = row_vec
         .par_iter()
         .zip(&col_vec)
@@ -421,6 +400,8 @@ pub(crate) fn arithmetize_matrix<F: PrimeField>(
         })
         .collect();
 
+    end_timer!(lde_evals_time);
+
     let interpolate_time = start_timer!(|| "Interpolating on K and B");
     let row_evals_on_K = EvaluationsOnDomain::from_vec_and_domain(row_vec, domain_k.clone());
     let col_evals_on_K = EvaluationsOnDomain::from_vec_and_domain(col_vec, domain_k.clone());
@@ -490,23 +471,26 @@ fn is_in_ascending_order<T: Ord>(x_s: &[T], is_less_than: impl Fn(&T, &T) -> boo
 
 /// This function converts a R1CS matrix from ginger-lib into the sparse matrix representation
 /// `Matrix` as used in this crate.
+/// The columns of the matrix are re-arranged to be coherent with the treatment of input variables
+/// as a subdomain of the full variable vector.
 fn to_matrix_helper<F: PrimeField>(
     matrix: &[Vec<(F, VarIndex)>],
-    num_input_variables: usize,
+    domain_h: &Box<dyn EvaluationDomain<F>>,
+    domain_x: &Box<dyn EvaluationDomain<F>>,
 ) -> SparseMatrix<F> {
     let mut new_matrix = Vec::with_capacity(matrix.len());
-    let domain_x = get_best_evaluation_domain::<F>(num_input_variables).unwrap();
     let domain_x_size = domain_x.size();
     for row in matrix {
         let mut new_row = Vec::with_capacity(row.len());
         for (fe, column) in row {
             let column = match column {
-                // public inputs correspond to the first columns
                 VarIndex::Input(i) => *i,
-                // private witnesses start right after
                 VarIndex::Aux(i) => domain_x_size + i,
             };
-            new_row.push((*fe, column))
+            let index = domain_h
+                .reindex_by_subdomain(domain_x_size, column)
+                .unwrap();
+            new_row.push((*fe, index))
         }
         new_matrix.push(new_row)
     }
@@ -537,11 +521,12 @@ fn balance_matrices<F: Field>(
 
 pub(crate) fn post_process_matrices<F: PrimeField>(
     cs: &mut ConstraintSystem<F>,
+    domain_h: &Box<dyn EvaluationDomain<F>>,
+    domain_x: &Box<dyn EvaluationDomain<F>>,
 ) -> Option<(SparseMatrix<F>, SparseMatrix<F>, SparseMatrix<F>)> {
-    balance_matrices(&mut cs.at, &mut cs.bt);
-    let a = to_matrix_helper(&cs.at, cs.num_inputs);
-    let b = to_matrix_helper(&cs.bt, cs.num_inputs);
-    let c = to_matrix_helper(&cs.ct, cs.num_inputs);
+    let a = to_matrix_helper(&cs.at, domain_h, domain_x);
+    let b = to_matrix_helper(&cs.bt, domain_h, domain_x);
+    let c = to_matrix_helper(&cs.ct, domain_h, domain_x);
     Some((a, b, c))
 }
 
@@ -611,16 +596,27 @@ mod test {
             let mut matrix =
                 build_random_matrix(num_inputs, num_aux, num_constraints, fill_factor, rng);
             let num_non_zero = matrix.iter().map(|row| row.len()).sum();
-            let (domain_h, domain_k, domain_x, domain_b) =
+            let (domain_h, domain_k, _domain_x, domain_b) =
                 IOP::build_domains(num_aux, num_inputs, num_constraints, num_non_zero).unwrap();
 
             let arithmetization =
-                arithmetize_matrix("m", &mut matrix, &domain_k, &domain_h, &domain_x, &domain_b)
-                    .unwrap();
+                arithmetize_matrix("m", &mut matrix, &domain_k, &domain_h, &domain_b).unwrap();
 
             // Random points for checking that dense and sparse representations of the matrix coincide
-            let x = F::rand(rng);
-            let y = F::rand(rng);
+            let x = loop {
+                let x = F::rand(rng);
+                // check that x does not belong to domain_k
+                if !domain_k.evaluate_vanishing_polynomial(x).is_zero() {
+                    break x;
+                }
+            };
+            let y = loop {
+                let y = F::rand(rng);
+                // check that y does not belong to domain_k
+                if !domain_k.evaluate_vanishing_polynomial(y).is_zero() {
+                    break y;
+                }
+            };
 
             // Evaluate M(x,y) using dense representation
             //   M(X,Y) = \sum_{i,j = 1...n} M_{i,j} * L(X, z_i) * L(Y, z_j)
@@ -634,9 +630,7 @@ mod test {
 
             for (i, row) in matrix.iter().enumerate() {
                 for &(val, j) in row {
-                    result_dense += val
-                        * lagrange_h_x[i]
-                        * lagrange_h_y[domain_h.reindex_by_subdomain(domain_x.size(), j).unwrap()];
+                    result_dense += val * lagrange_h_x[i] * lagrange_h_y[j];
                 }
             }