CKKS同态加密通用函数近似方法和openFHE实现

摘要

同态加密可以直接在密文上进行运算，尤其是CKKS，可以直接在实数的密文上进行运算。服务器可以利用强大的计算能力，在不泄露用户隐私的情况下，为用户提供便捷的外包运算服务。然而，CKKS只能进行算术运算（多项式函数），无法直接计算很多复杂的函数。本文介绍了如何使用CKKS计算近似任意函数，并介绍了利用openFHE进行代码实现。

多项式近似

根据魏尔斯特拉斯近似定理，定义在闭区间上的连续函数可被多项式函数任意接近地一致近似。也就是当多项式的次数足够高的时候，在闭区间上的任意一点，我们都可以得到非常好的近似结果。

因此，使用同态加密解决问题的一般流程如下：

多项式近似在同态加密的应用中具有非常重要的地位。

常用的多项式近似方法有：Remes近似、切比雪夫多项式近似、泰勒展开、最小二乘法等。其中切比雪夫多项式逼近可以最大限度地降低龙格现象，也就是不会在某一点上的误差特别大。

当然，并不是所有函数都可以这样做，因为很多函数不是连续的，比如符号函数，取模函数等。这些函数在某些区间上可以得到很好的近似，但是在断点附近，会产生很大的误差。这种函数通常的做法是，使用复合函数（逼近符号函数），或者首先使用三角函数等基本的连续函数去逼近（如取模函数），然后再使用多项式逼近。

openFHE实现

openFHE实现了一个使用切比雪夫多项式逼近的函数接口，官方例子。

密文上下文对象的 EvalChebyshevFunction 函数实现了任意函数的切比雪夫多项式近似。其输入为一个函数，需要计算的密文，输入区间的下界和上界，切比雪夫多项式的次数。

作为输入的函数，其参数是一个 double 类型的变量，输出是 double ，下面是一个例子：

double func(double x ){
    return std::sin(x);
}

输入密文需要留有足够的乘法深度。乘法的深度和切比雪夫多项式的次数有关系，下面是openFHE官方给的表格：

具体的多项式计算采用了树结构，所以多项式度数和乘法深度大约是对数关系。

近似是否准确，主要取决于两个因素，一个是输入的区间范围，一个是多项式的次数。多项式的次数越高，那么计算需要的时间越长。理论上的近似精度会更高，但这只是从全局的误差来看，具体的某个点的误差则不确定。

输入区间的范围对近似精度的影响非常大，所以，在选择区间的时候，需要格外注意。

下面是我从openFHE的GitHub上下载的代码，小改后的版本：

#include<iostream>
#include"openfhe.h"
//The functions or classes of OpenFHE are in the namespace lbcrypto
using namespace lbcrypto;
using namespace std;

double modfunc(double x){
    return sqrt(x);
}

void EvalFunctionExample() {
    std::cout << "--------------------------------- EVAL SQUARE ROOT FUNCTION ---------------------------------"
              << std::endl;
    CCParams<CryptoContextCKKSRNS> parameters;

    // We set a smaller ring dimension to improve performance for this example.
    // In production environments, the security level should be set to
    // HEStd_128_classic, HEStd_192_classic, or HEStd_256_classic for 128-bit, 192-bit,
    // or 256-bit security, respectively.
    parameters.SetSecurityLevel(HEStd_128_classic);
    //parameters.SetRingDim(1 << 10);
#if NATIVEINT == 128
    usint scalingModSize = 78;
    usint firstModSize   = 89;
#else
    usint scalingModSize = 50;
    usint firstModSize   = 60;
#endif
    parameters.SetScalingModSize(scalingModSize);
    parameters.SetFirstModSize(firstModSize);

    // Choosing a higher degree yields better precision, but a longer runtime.
    uint32_t polyDegree=50;
    // The multiplicative depth depends on the polynomial degree.
    // See the FUNCTION_EVALUATION.md file for a table mapping polynomial degrees to multiplicative depths.
    uint32_t multDepth = 13;

    parameters.SetMultiplicativeDepth(multDepth);
    CryptoContext<DCRTPoly> cc = GenCryptoContext(parameters);
    cc->Enable(PKE);
    cc->Enable(KEYSWITCH);
    cc->Enable(LEVELEDSHE);
    // We need to enable Advanced SHE to use the Chebyshev approximation.
    cc->Enable(ADVANCEDSHE);

    auto keyPair = cc->KeyGen();
    // We need to generate mult keys to run Chebyshev approximations.
    cc->EvalMultKeyGen(keyPair.secretKey);

    std::vector<std::complex<double>> input{1, 2, 3, 4, 5, 6, 7, 8, 9};
    size_t encodedLength = input.size();
    Plaintext plaintext  = cc->MakeCKKSPackedPlaintext(input);
    auto ciphertext      = cc->Encrypt(keyPair.publicKey, plaintext);

    double lowerBound = 0;
    double upperBound = 10;

    // We can input any lambda function, which inputs a double and returns a double.
    //auto result = cc->EvalChebyshevFunction([](double x) -> double { return std::sqrt(x); }, ciphertext, lowerBound,upperBound, polyDegree);
    auto result = cc->EvalChebyshevFunction(modfunc, ciphertext, lowerBound,
                                            upperBound, polyDegree);
    Plaintext plaintextDec;
    cc->Decrypt(keyPair.secretKey, result, &plaintextDec);
    plaintextDec->SetLength(encodedLength);

    std::vector<std::complex<double>> expectedOutput(
        {sqrt(1), sqrt(2), sqrt(3), sqrt(4), sqrt(5), sqrt(6), sqrt(7), sqrt(8), sqrt(9)});
    std::cout << "Expected output\n\t" << expectedOutput << std::endl;

    std::vector<std::complex<double>> finalResult = plaintextDec->GetCKKSPackedValue();
    std::cout << "Actual output\n\t" << finalResult << std::endl << std::endl;
}

int main(){
    EvalFunctionExample();
}

下面是lowerBound=0, upperBound=10时的结果：

下面是lowerBound=0, upperBound=20时的结果：

下面是lowerBound=0, upperBound=100时的结果：