2021SC@SDUSC

看一看估计标准偏差过程是怎么样的

// Estimate standard deviation using the imaginary part of decoded vector z

// Compute m(X) - m(1/X) as a proxy for z - Conj(z) = 2*Im(z)

// vec is m(X) corresponding to z

// conjugate is m(1/X) corresponding to Conj(z)

首先可以关注一下这段注释

M(x) - M(1/x) as a Proxy..

 using the imaginary part of decoded vector z

z - Conj(z) = 2*Im(z)

也就是说需要利用共轭的部分 其实这也就是共轭的定义式了

z - Conj(z) = 2*Im(z)

但是具体的含义不太清楚

可以看一下以下这段代码

不知道cck的编码方式是否要先虚数话然后再行具体的编码

#include "encoding/ckkspackedencoding.h"

#include "math/dftransfrm.h"

namespace lbcrypto {

std::vector<std::complex<double>> Conjugate(

const std::vector<std::complex<double>> &vec) {

uint32_t n = vec.size();

std::vector<std::complex<double>> result(n);

for (size_t i = 1; i < n; i++) {

result[i] = {-vec[n - i].imag(), -vec[n - i].real()};

}

result[0] = {vec[0].real(), -vec[0].imag()};

return result;

}

// Estimate standard deviation using the imaginary part of decoded vector z

// Compute m(X) - m(1/X) as a proxy for z - Conj(z) = 2*Im(z)

// vec is m(X) corresponding to z

// conjugate is m(1/X) corresponding to Conj(z)

double StdDev(const std::vector<std::complex<double>> &vec,

const std::vector<std::complex<double>> &conjugate) {

uint32_t Nh = vec.size();

  // ring dimension

uint32_t n = Nh * 2;

  // extract the complex part using identity z - Conj(z) == 2*Im(z)

  // here we actually compute m(X) - m(1/X) corresponding to 2*Im(z).

  // we only need first Nh/2 + 1 components of the imaginary part

  // as the remaining Nh/2 - 1 components have a symmetry

  // w.r.t. components from 1 to Nh/2 - 1

std::vector<std::complex<double>> complexValues(Nh / 2 + 1);

for (size_t i = 0; i < Nh / 2 + 1; ++i) {

complexValues[i] = vec[i] - conjugate[i];

}

  // Calculate the mean

auto mean_func = [](double accumulator, const std::complex<double> &val) {

return accumulator + (val.real() + val.imag());

};

  // use the symmetry condition

double mean =

2 * std::accumulate(complexValues.begin() + 1,

complexValues.begin() + Nh / 2, 0.0, mean_func);

  // and then add values at indices 0 and Nh/2

mean += complexValues[0].imag();

mean += 2 * complexValues[Nh / 2].real();

  // exclude the real part at index 0 as it is always 0

mean /= static_cast<double>(n) - 1.0;

  // Now calculate the variance

auto variance_func = [&mean](double accumulator,

const std::complex<double> &val) {

return accumulator + (val.real() - mean) * (val.real() - mean) +

(val.imag() - mean) * (val.imag() - mean);

};

  // use the symmetry condition

double variance =

2 * accumulate(complexValues.begin() + 1, complexValues.begin() + Nh / 2,

0.0, variance_func);

  // and then add values at indices 0 and Nh/2

variance +=

(complexValues[0].imag() - mean) * (complexValues[0].imag() - mean);

variance += 2 * (complexValues[Nh / 2].real() - mean) *

(complexValues[Nh / 2].real() - mean);

  // exclude the real part at index 0 as it is always 0

variance /= static_cast<double>(n) - 2.0;

  // scale down by 2 as we have worked with 2*Im(z) up to this point

double stddev = 0.5 * std::sqrt(variance);

return stddev;

}

std::vector<DCRTPoly::Integer> CKKSPackedEncoding::CRTMult(

const std::vector<DCRTPoly::Integer> &a,

const std::vector<DCRTPoly::Integer> &b,

const std::vector<DCRTPoly::Integer> &mods) {

std::vector<DCRTPoly::Integer> result(mods.size());

for (usint i = 0; i < a.size(); i++) {

result[i] = a[i].ModMulFast(b[i], mods[i]);

}

return result;

}

但是由以上可见 这个实现过程并不容易 (不过这里有一个问题 那么是否只有cck才使用这样的编码方式呢 如果不是这样的话共轭计算可以作为core/lib/encoding下的一个公共的库函数的)

再说回这个具体的过程 最终的成果是一个实数 而输入则是一个向量

估计标准误差（Se）是说明实际值与其估计值之间相对偏离程度的指标，主要用来衡量 回归方程 的代表性。 估计标准误差的值越小，则 估计量 与其真实值的近似误差越小，但不能认为估计量与真实值之间的 绝对误差 就是估计标准误差。

这是一个统计学中的计算概念 而这里实现了它

再来看看接下来的部分

palisade作为一个跨平台的库 这里也进行了基于宏的机型适配

不过这两段方法之间会不会由相当可以提取重用逻辑的地方呢

还是由于本地Int长度的差别其中的优化和计算方式差别较大呢

#if NATIVEINT == 128

bool CKKSPackedEncoding::Encode() {

#else  // NATIVEINT == 64

bool CKKSPackedEncoding::Encode()

那么ccks具体的编码方式是什么样的呢

首先先确定该条是否已经进行了编码

然后确认生一个多项式环

Nh已经出现了多次 但是目前尚未知道其具体含义

被加密的对象只能是DCRTPoly

一种特殊的多项式

这的概念有两篇很好的blog已经做了详尽的分析

(32条消息) 2021SC@SDUSC PALISADE开源库（二）CKKS讲解系列（一）普通编码和解码_小河豚oO的博客-CSDN博客

(32条消息) 2021SC@SDUSC PALISADE开源库（二）CKKS讲解系列（二）完整编码和解码_小河豚oO的博客-CSDN博客_palisade库

if (this->isEncoded) return true;

uint32_t ringDim = GetElementRingDimension();

uint32_t Nh = (ringDim >> 1);

std::vector<std::complex<double>> inverse = this->GetCKKSPackedValue();

  // clears all imaginary values as CKKS for complex numbers

for (size_t i = 0; i < inverse.size(); i++) inverse[i].imag(0.0);

inverse.resize(Nh);

if (this->typeFlag == IsDCRTPoly) {

DiscreteFourierTransform::FFTSpecialInv(inverse);

uint64_t pBits = encodingParams->GetPlaintextModulus();

uint32_t precision = 52;

double powP = std::pow(2, precision);

int32_t pCurrent = pBits - precision;

    // the idea is to break down real and imaginary parts

    // expressed as input_mantissa * 2^input_exponent

    // into (input_mantissa * 2^52) * 2^(p - 52 + input_exponent)

    // to preserve 52-bit precision of doubles

    // when converting to 128-bit numbers

std::vector<__int128> temp(2 * Nh);

for (size_t i = 0; i < Nh; ++i) {

      // Check for possible overflow in llround function

int32_t n1 = 0;

      // extract the mantissa of real part and multiply it by 2^52

double dre =

static_cast<double>(std::frexp(inverse[i].real(), &n1) * powP);

int32_t n2 = 0;

      // extract the mantissa of imaginary part and multiply it by 2^52

double dim =

static_cast<double>(std::frexp(inverse[i].imag(), &n2) * powP);

if (is128BitOverflow(dre) || is128BitOverflow(dim)) {

PALISADE_THROW(math_error, "Overflow, try to decrease scaling factor");

}

int64_t re64 = std::llround(dre);

int32_t pRemaining = pCurrent + n1;

__int128 re = 0;

if (pRemaining < 0) {

re = re64 >> (-pRemaining);

} else {

__int128 pPowRemaining = ((__int128)1) << pRemaining;

re = pPowRemaining * re64;

}

int64_t im64 = std::llround(dim);

pRemaining = pCurrent + n2;

__int128 im = 0;

if (pRemaining < 0) {

im = im64 >> (-pRemaining);

} else {

__int128 pPowRemaining = ((int64_t)1) << pRemaining;

im = pPowRemaining * im64;

}

temp[i] = (re < 0) ? Max128BitValue() + re : re;

temp[i + Nh] = (im < 0) ? Max128BitValue() + im : im;

if (is128BitOverflow(temp[i]) || is128BitOverflow(temp[i + Nh])) {

PALISADE_THROW(math_error, "Overflow, try to decrease scaling factor");

}

}

const shared_ptr<ILDCRTParams<BigInteger>> params =

this->encodedVectorDCRT.GetParams();

const std::vector<std::shared_ptr<ILNativeParams>> &nativeParams =

params->GetParams();

for (size_t i = 0; i < nativeParams.size(); i++) {

NativeVector nativeVec(ringDim, nativeParams[i]->GetModulus());

FitToNativeVector(temp, Max128BitValue(), &nativeVec);

NativePoly element = this->GetElement<DCRTPoly>().GetElementAtIndex(i);

element.SetValues(

nativeVec, Format::COEFFICIENT);  // output was in coefficient format

this->encodedVectorDCRT.SetElementAtIndex(i, element);

}

usint numTowers = nativeParams.size();

std::vector<DCRTPoly::Integer> moduli(numTowers);

for (usint i = 0; i < numTowers; i++) {

moduli[i] = nativeParams[i]->GetModulus();

}

DCRTPoly::Integer intPowP = NativeInteger(1) << pBits;

std::vector<DCRTPoly::Integer> crtPowP(numTowers, intPowP);

auto currPowP = crtPowP;

    // We want to scale temp by 2^(pd), and the loop starts from j=2

    // because temp is already scaled by 2^p in the re/im loop above,

    // and currPowP already is 2^p.

for (size_t i = 2; i < depth; i++) {

currPowP = CKKSPackedEncoding::CRTMult(currPowP, crtPowP, moduli);

}

if (depth > 1) {

this->encodedVectorDCRT = this->encodedVectorDCRT.Times(currPowP);

}

this->GetElement<DCRTPoly>().SetFormat(Format::EVALUATION);

scalingFactor = pow(scalingFactor, depth);

} else {

PALISADE_THROW(config_error, "Only DCRTPoly is supported for CKKS.");

}

this->isEncoded = true;

return true;

那我们具体再来看看这64位的和128位的差别在于哪里

if (is64BitOverflow(dre) || is64BitOverflow(dim)) {

        // IFFT formula:

        // x[n] = (1/N) * \Sum^(N-1)_(k=0) X[k] * exp( j*2*pi*n*k/N )

        // n is i

        // k is idx below

        // N is inverse.size()

        //

        // In the following, we switch to original data domain,

        // and we identify the component that has the maximum

        // contribution to the values in the iFFT domain. We do

        // this to report it to the user, so they can identify

        // large inputs.

由此可以看到 大概的差别是在计算方式上 也就是说 对于大整数的计算上和精度较高的计算上

需要进行一个映射以适应精度要求