[email protected]
看一看估计标准偏差过程是怎么样的
// 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条消息) [email protected] PALISADE开源库（二）CKKS讲解系列（一）普通编码和解码_小河豚oO的博客-CSDN博客
(32条消息) [email protected] 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.
由此可以看到 大概的差别是在计算方式上 也就是说 对于大整数的计算上和精度较高的计算上
需要进行一个映射以适应精度要求
![C语言第四章自述2第四章 选择结构程序设计[图]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)