Polynomial Evaluation on BGV

[Th0masAndy]

func ObliviousPolynomialEvaluationBGV(params bgv.Parameters, degree int, ct *rlwe.Ciphertext, pre_power *rlwe.Ciphertext, coeffs_pt []*rlwe.Plaintext, evaluators []*bgv.Evaluator, encryptors []*rlwe.Encryptor, decryptor *rlwe.Decryptor, encoder *bgv.Encoder, targetLevel int) *rlwe.Ciphertext {
	var err error
	square_root := int(math.Ceil(math.Sqrt(float64(degree))))
	zero_cnt := 0
	blocks := make([]*rlwe.Ciphertext, square_root)
	for i := range blocks {
		blocks[i] = getItem(zero_pool)
		zero_cnt++
	}
	logger.Printf("=========== Evaluate Poly of degree %d ===========\n\n", degree)
	start := time.Now()
	small_base_power_num := int(math.Ceil(math.Log2(float64(square_root))))
	small_base_power := make([]*rlwe.Ciphertext, small_base_power_num+2)
	small_base_power[0] = getItem(one_pool) // one
	small_base_power[1] = ct                // x**1
	for i := 2; i < small_base_power_num+2; i++ {
		tmp, err := evaluators[0].MulRelinNew(small_base_power[i-1], small_base_power[i-1])
		if err != nil {
			panic(err)
		}
		err = evaluators[0].Rescale(tmp, tmp)
		if err != nil {
			panic(err)
		}
		small_base_power[i] = tmp
	}
	// power of x 0 1 2 4 8
	small_power_num := square_root
	small_power := make([]*rlwe.Ciphertext, small_power_num+1)
	small_power[0] = small_base_power[0] //one
	small_power[1] = small_base_power[1] //x**1
	wg.Add(small_power_num - 1)
	for i := 2; i < small_power_num+1; i++ {
		go func(i int) {
		index := i
		var one *rlwe.Ciphertext
		for j := 0; j <= small_base_power_num; j++ {
			if index&1 == 1 {
				if one == nil {
					one = small_base_power[j+1]
				} else {
					one, err = evaluators[i-2].MulRelinNew(small_base_power[j+1], one)
					if err != nil {
						panic(err)
					}
					err := evaluators[i-2].Rescale(one, one)
					if err != nil {
						panic(err)
					}
				}
			}
			index = index >> 1
		}
		small_power[i] = one
		res_mat := decryptor.DecryptNew(small_power[i])
		res := make([]uint64, params.MaxSlots())
		encoder.Decode(res_mat, res)
		fmt.Println("small pow", i, "level", small_power[i].Level(), res[:10])
		 wg.Done()
		 }(i)
	}
	wg.Wait()

	big_base_power_num := int(math.Ceil(math.Log2(math.Ceil(math.Sqrt(float64(degree)) - 1))))
	big_base_power := make([]*rlwe.Ciphertext, big_base_power_num+2)
	big_base_power[0] = small_base_power[0] //x**0
	if pre_power != nil {
		big_base_power[1] = pre_power //pre compute power x**16
	} else {
		big_base_power[1] = small_power[small_power_num] //the last power
	}
	for i := 2; i < big_base_power_num+2; i++ {
		tmp, err := evaluators[0].MulRelinNew(big_base_power[i-1], big_base_power[i-1])
		if err != nil {
			panic(err)
		}
		err = evaluators[0].Rescale(tmp, tmp)
		if err != nil {
			panic(err)
		}
		big_base_power[i] = tmp
	}
	// power of x**16 0 1 2 4 8 16

	big_power_num := square_root
	big_power := make([]*rlwe.Ciphertext, big_power_num)
	big_power[0] = big_base_power[0] //x**0
	big_power[1] = big_base_power[1] //x**16

	wg.Add(big_power_num - 2)

	for i := 2; i < big_power_num; i++ {
		go func(i int) {
		index := i
		var one *rlwe.Ciphertext
		for j := 0; j <= big_base_power_num; j++ {
			if index&1 == 1 {
				if one == nil {
					one = big_base_power[j+1]
				} else {
					one, err = evaluators[i-2].MulRelinNew(big_base_power[j+1], one)
					if err != nil {
						panic(err)
					}
					err := evaluators[i-2].Rescale(one, one)
					if err != nil {
						panic(err)
					}
				}
			}
			index = index >> 1
		}
		big_power[i] = one
		res_mat := decryptor.DecryptNew(big_power[i])
		res := make([]uint64, params.MaxSlots())
		encoder.Decode(res_mat, res)
		fmt.Println("big pow", i*square_root, "level", big_power[i].Level(), res[:10])
		 wg.Done()
		 }(i)
	}
	 wg.Wait()

	// compute every blocks
	coeffs := make([][]uint64, degree)
	for i := 0; i < degree; i++ {
		t := make([]uint64, params.N())
		for j := range t {
			t[j] = uint64(i + 1)
		}
		coeffs[i] = t
	}
	 wg.Add(big_power_num)
	for i := 0; i < big_power_num; i++ {
		 go func(i int) {
		level := targetLevel + 1
		sub_res_tmp := rlwe.NewCiphertext(params, 1, level)
		sub_res_tmp.Scale = params.DefaultScale().Div(big_power[i].Scale).Mul(rlwe.NewScale(params.RingQ().AtLevel(level).SubRings[level].Modulus))
		zero_cnt++
		for j := 0; j < small_power_num; j++ {
			if i*small_power_num+j == degree {
				break
			}
			evaluators[i].MulThenAdd(small_power[j], coeffs[i*small_power_num+j], sub_res_tmp)
		}
		fmt.Println(sub_res_tmp.Level(), big_power[i].Level())
		blocks[i], err = evaluators[i].MulRelinNew(sub_res_tmp, big_power[i])
		if err != nil {
			panic(err)
		}
		evaluators[i].Rescale(blocks[i], blocks[i])

		fmt.Println("level:", blocks[i].Level(), "block scale:", blocks[i].Scale.Uint64())

		 wg.Done()
		 }(i)
	}
	wg.Wait()

	res := blocks[0]

	for i := range blocks[1:] {
		res, err = evaluators[0].AddNew(res, blocks[i])

		if err != nil {
			panic(err)
		}
		res_mat := decryptor.DecryptNew(res)
		res_pt := make([]uint64, params.MaxSlots())
		encoder.Decode(res_mat, res_pt)
		fmt.Println("blocks", i, res_pt[:10], res.Scale.Uint64())
	}

	dt := time.Since(start)
	logger.Printf("> Eval Poly time %v \n\n", dt)
	logger.Printf("> Zero pop count %v \n\n", zero_cnt)

	fmt.Println("Output Scale:", res.Scale.Uint64())
	return res
}

I implemented the Paterson Stockmeyer algorithm and get the right result.
For example, I want to compute a polynomial f(x) of degree 255, fisrt compute small steps $x^0,x^1.....x^{15}$ and big steps $x^{16},x^{32},x^{48}..... x^{240}$ and then compute the blocks.

The block 0 is $a_0 x^0+a_1 x^1.....a_{15} x^{15}$ , here $a_i$ is the coeffs
the block 1 is $x^{16} (a_{16} + a_{17} x^1 ..... + a_{31}x^{15})$ and $block_i$ is the same
and finally I add these blocks together to get $f(x) = block_0 + block_1 .... + block_{15}$ .
The result of f(x) is right and it costs 9 multiplications to compute f(x).
I choose param

PQ := bgv.ParametersLiteral{
		LogN:             15,
		LogQ:             []int{40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40},
		LogP:             []int{40, 40, 40},
		PlaintextModulus: 65537,
	}

and I can compute 18 multiplications at most for this param.
I'm supposed to be able to perform 9 multiplications based on f(x) but I found out I can only do 7 multiplications.

Will my algorithm significantly reduce the noise budget?

[Pro7ech]

Hi @Th0masAndy, thank you for your question and interest in the Lattigo library. I have converted your issue to a discussion, as I believe this is more appropriate (see the README for what issue should be used for).

The library already provides a polynomial evaluation algorithm based on the Paterson Stockmeyer with power of two decomposition. So I believe that your re-implementation is to add concurrency, is this correct?

[Pro7ech]

The multiplicative depth of the algorithm implemented in Lattigo is already optimal : $\lceil\log_2(degree+1)\rceil$, it is proven in Section 3.3 of https://eprint.iacr.org/2020/1203.

Edit: ok I see what you meant, having the client pre-compute the power basis. This also possible and already provided by Lattigo by using the PowerBasis struct, and pre-compute the necessary powers and giving it as input to the algorithm instead of the original ciphertext.

Then what remains to do is to add concurrency. For this I suggest that you start from this branch: https://github.com/tuneinsight/lattigo/blob/dev_release_v4.2.0/circuits/polynomial_evaluator.go which is our next release and provides an inlined version of the algorithm, with interfaces, which can be easily parallelized. It is not stable yet, but should provide you with a good starting point.

[AsterNighT]

Thank you for the advice. We would look in to it. But our first target in raising this discussion is to seek insight on the reason why our implementation is inefficient and sometimes even incorrect. I apologize for not making the point clearer at the very beginning.

I would like to provide some extra clue discovered in our debugging, the code does runs on the dev_release_v4.2.0 branch. We believe it has something to do with the scale of the ciphertexts:

• Evaluator.Add(c1,c0,c0) and Evaluator.Add(c0,c1,c0) has very different result. The former is incorrect, while the latter is correct. They does not seem to run out of noise budget.
• Multiplying the sum of a serial of ciphertexts with different level and scale often yields incorrect results. But by aligning the scale manually, the results are more likely to be correct.

[Pro7ech]

Ok, it's clear now. The thing with BGV is that scale alignment, if not done carefully, can consumes up to a full level of noise budget. And the reason why the scale changes depending on the order of the inputs is because the scale matching factors are computed as (r0, r1) = f(input_1, input_2), so switching the order of the input will switch the order of r0, r1, and the result will be scaled by r0 or r1 depending on the order of the inputs.

[AsterNighT]

Thank you, that's most informative for us. Still I'm not quite sure what is the exact reason scaling would consume noise budget. I once read about the paper Fully Homomorphic Encryption without Bootstrapping but I still didn't quite get the idea behind all these scaling issues. Is there something I missed in the paper or the scaling is just a more intractable issue in implementation?

Besides, I wonder if you could provide us(and those who are not so familiar with bgv scale management) some advice on how to avoid such noise budget loss? I notice in the dev-poly-eval branch there is a dummy evaluator for precalculating scales. Is it the common practice when handling BGV ciphertext?

[Pro7ech]

Scaling can consume up to one level because you need to multiply one of the ciphertext by the difference. However, you can minimize that by doing the following. Given two ciphertexts at scale $c_0$ and $c_1$ with $c_0 \neq c_1$, you can find $r_0$ and $r_1$ such that $c_0r_0 = c_1r_1$ and $|r_0+r_1|$ is minimal. This can be done using the extended Euclidean algorithm (https://github.com/tuneinsight/lattigo/blob/98b412dc0e50068352e9b82b5890896d4b2a5af0/bgv/evaluator.go#L1412)

The polynomial evaluation implemented in Lattigo computes the exact scaling factors for each of the polynomial coefficients such that no scaling operation is needed when adding two ciphertexts and such that the output ciphertext has the desired scale.

This is done by first simulating the evaluation and computing what scale each ciphertext should have at each step of the evaluation, propagating theses scales down the evaluation tree, until we reach the plaintext coefficients. These plaintext coefficients are then scaled by the appropriate amount. The resulting behavior is that all intermediate additions are done between ciphertexts that have the exact scale. And the final ciphertext has the desired scale.

The original idea behind this algorithm comes from Section 3 of this paper: https://eprint.iacr.org/2020/1203