Privacy Preserving Recursive Least Squares

107 days ago by ktjell12@student.aau.dk

import numpy as np def o(x): return np.array(x, dtype = 'object') #Creates shares of secrets using Shamir's secret sharing scheme. def secretsharing(x, t, n): shares = [] c = [F.random_element() for i in range(t)] for i in range(1, n+1): s = int(x) for j in range(1,t+1): s = s + o(c[j-1]) * o(pow(i,j)) shares.append(s) return o(shares) #Creates the "recombination"-vector used to reconstruct a secret from its shares. def basispoly(n): r = [] C = range(1,n+1) for i in range(1,n+1): c = [k for k in C if k != i] p = 1 for j in range(n-1): p *= -F(c[j]) / (F(i)-F(c[j])) r.append(F(p)); return r # reconstruct secret. def rec(x): res = 0 n = len(x) y = basispoly(n) for i in range(len(x)): res += o(x[i]) * o(y[i]) return res 
       
#Create a Beaver's triplet. def triplet(n, t): r = basispoly(n) matrixA = np.zeros((n,n), dtype = 'object') matrixB = np.zeros((n,n), dtype = 'object') for i in range(n): a = F.random_element() b = F.random_element() matrixA[:,i] = secretsharing( a, t, n) matrixB[:,i] = secretsharing(b, t, n) sharesA = np.sum(matrixA, 1) sharesA = [F(i) for i in sharesA] sharesB = np.sum(matrixB, 1) sharesB = [F(i) for i in sharesB] sharesC = mult(sharesA,sharesB) #[F(i*j) for (i,j) in zip(sharesA, sharesB)] return [sharesA, sharesB, sharesC] #Secure mulitplication of two secrets. -Resulting polynomial is still degree t. def mult(a,b): r = basispoly(n) h = o(a)*o(b) hs = [] for i in range(n): hs.append(secretsharing(h[i],t,n)) s = 0 for i in range(n): s+= o(r[i]) * o(hs[i]) return o(s) # Return a random bit unknown to all parties. def randomBit(n, t): rsq = 0 while rsq == 0: #If r**2 = 0, the parties retry r_shares = triplet(n,t)[0] #The parties compute an unknown random number, r, and computes r**2 rsq_shares = mult(r_shares,r_shares) rsq = rec(rsq_shares) #r**2 is revealed r_prime = sqrt(rsq) return (F(1)/2 * (o(r_prime**(-1)) * o(r_shares) + 1)) #Return both sharings of a random number and the bitwise sharing of the same random number. def randomBitNumber(n, t, l): sharings = [] for i in range(l): sharings.append(randomBit(n,t)) r_shares = 0 for j in range(l): r_shares += o(2**j) * o(sharings[j]) return sharings, r_shares 
       
##################################################################################################### # This cell implementes SMPC functions. # # For instance: # # - Multiplication # # - Inverse field element # # - Bitwise operations (OR, AND, XOR) # # - Bit decomposition # # - Comparison (bitwise Less-Than) # ##################################################################################################### #Secure mulitplication of two secrets. -Resulting polynomial is still degree t. def mult(a,b): r = basispoly(n) h = o(a)*o(b) hs = [] for i in range(n): hs.append(secretsharing(h[i],t,n)) s = secretsharing(0,t,n) for i in range(n): s+= o(r[i]) * o(hs[i]) return np.array(s,dtype = 'object') #Computing the inverse of a field element. def inverse(x, r): w_s = mult(x,r) w = rec(w_s) return w**(-1) * o(r) #Bitwise AND operation def AND(a): p = a[0] for i in range(1,len(a)): p = mult(p, a[i]) return p #Bitwise XOR operation def xor(a1,b1): d = o(a1) - o(b1) return mult(d,d) #Bitwise OR operation def OR(a): p = 1 - a[0] for i in range(1,len(a)): p = mult(p, (1 - a[i])) return 1 - p #Converting bitwise sharing's of a secret to ordinary sharings. def dec(a): s = 0 for i in range(len(a)): s+= o(a[i]) * o(2**i) return s #Computing the truth of whether a < b. def lessthan(a,b): e = [] for i in range(len(a)): e.append(xor(a[i],b[i])) f = [] for i in range(len(e)-1): f.append(OR(e[i:])) f.append(e[-1]) g = [] for i in range(0,len(f)-1): g.append(f[i] - f[i+1]) g.append(f[-1]) h = [] for i in range(len(g)): h.append(mult(g[i], b[i])) y = 0 for j in range(len(h)): y += h[j] return y #Computation of BitLen without bitwise decomposition. #Computes the bit length of a secret def bitLen(F, d, r, rb, n, t, l): ### ADD RANDOMNESS rb.append(secretsharing(0,t,n)) es = o(d) + o(r) e = rec(es) ebit = list(bin(e)[2:].rjust(l+1,'0')) ebit = ebit[::-1] ebs = [secretsharing(b,t,n) for b in ebit] ### Create H ep = [xor(i,j) for i,j in zip(ebs,rb)] E = [] for i in range(len(ep)): E.append(OR(ep[i:])) E[-1] = ep[-1] h = [i - mult(i,j) for i,j in zip (rb, ebs)] hp = [] for i in range(len(h)): if i == len(h)-1: Es = 0 else: Es = E[i+1] hp.append(h[i] + 1 - Es) H = [] for i in range(len(hp)): H.append(AND(hp[i:])) ### Calc Hld Hp = [] for i in range(len(H)): Hp.append(xor(H[i],1)) e0 = [] r0 = [] for i in range(len(Hp)): e0.append(mult(Hp[i],ebs[i])) r0.append(mult(Hp[i], rb[i])) B1 = lessthan(e0,r0) H_aim = [1-i for i in B1] ### Trans H into D D_ = [] H2 = H[:-1] H2.insert(0,secretsharing(0,t,n)) N = [1-i for i in H_aim] for i in range(len(H)): D_.append(mult(H_aim, H2[i]) + mult(N,H[i])) D = [1-i for i in D_] y = [] for i in range(len(D[0])): b = [] for j in range(len(D)): b.append(D[j][i]) y.append(dec(b)) return o(y) def binSub(a,b,l): y = [xor(a[0], b[0])] c = mult(1 - a[0], b[0]) for i in range(1,l): cp = mult(b[i], c) e = xor(b[i], c) y.append(xor(a[i], e)) c = cp + mult(1-a[i], e) return y def bitDec(a,l): rbits, r = randomBitNumber(n,t,l) rbits.append(secretsharing(0, t,n)) d = a + r dp = rec(d) db = bin(int(dp))[2:].rjust(l+1, '0') db = [int(i) for i in db] db = db[::-1] dpb = bin(int(dp) + F.cardinality())[2:].rjust(l+1, '0') dpb = [int(i) for i in dpb] dpb = dpb[::-1] b = lessthan(db, rbits) s = [b * dpb[i] + (1-b) * db[i] for i in range(l+1)] return binSub(s, rbits, l+1) #Divide secret with 2**k using bit decomposition (binary shift) def reScale(x, k, l): xb = bitDec(x, l) res = xb[k:] return dec(res) #Dividing a secret with a public integer def down(F, q, k,t,n): #while True: # r = F.random_element() # if r < 10*dot(F,basis,q): # break r = F(900)# r_T = np.trunc(np.array(r,dtype = float) *k**(-1)) r_s = secretsharing(r,t,n) r_Ts = secretsharing(F(r_T),t,n) z_s = q + r_s z = rec(z_s) z_T = np.trunc(np.array(z, dtype = float) * k**(-1)) z_Ts = secretsharing(F(z_T),t,n) return z_Ts - r_Ts #Integer division def integerDiv(F, d, nom, n, t, l, k): r1,r2,c = triplet(n, t) rr = 0 while rr == 0: rbits, r = randomBitNumber(n, t, l) rr = rec(r) D = bitLen(F, d, r, rbits, n, t, l) shift = 1 + o(D) #BitLen shiftInv = inverse(shift, r) a1 = o(shift) - o(d) p = mult(a1,shiftInv) prod = secretsharing(1,t,n) s = o(n*[0]) for i in range(1,l): prod = mult(prod, p) s = s + o(prod) s = s + 1 a_hat = 2**k * o(shiftInv) a_hat = mult(a_hat,s) q_hat = mult(nom,a_hat) q = reScale(q_hat, k, 240) #down(F, q_hat, 2**k,t,n) r = nom - mult(d,q) rd = bitDec(r+d, 2*l) ep = lessthan(rd, bitDec(d+d, 2*l)) ep = (ep - 1) * -1 em = lessthan(rd, bitDec(d, 2*l)) y = q + ep - em return y ###############################MAIN ########################################## ##### CONSTANTS ################################ set_random_seed(1) np.random.seed(1) m = random_prime(2**250, None, 2**247) # Prime number, large enough that the probability of chossing a certain number is # negeligable F = GF(m) # Finite field with m elements n = 4 # Number of parties t = 1 # Degree of polynomial l = 4 # bit length of some secrets 
       
########################################################################### #This cells defines functions that are more specifically used to # #compute the recursive least squares equations securely. # ########################################################################### #To be used when result is revealed. def less(a,b): return int( rec(a) < rec(b) ) # Recontruct a vector def reVector(b, n): Nvector = [] for i in range(2): liste = [] for j in range(n): liste.append(b[j][i][0]) Nvector.append(F(rec(liste)))#[b[0][1][i][0],b[1][1][i][0],b[2][1][i][0],b[3][1][i][0]])) return matrix(2,1,Nvector) #Integer-Least-Squares def ILS(F, P, x, w, y, C, t, n, i): g = matVec(P,x) e0 = vecVec(x,w) e0 = down(F,e0,C,t,n) e = [k-j for k,j in zip(y, e0)] a = [mult(k,e) for k in g] w = [i+j for i,j in zip(w,a)] return w #Estimate P-matrix def estP(F,P, x, C,t,n): xP = vecMat(x,P) a1 = vecVec(xP, x) denom = C + a1 Px = matVec(P,x) Part = VecVec(Px, xP) return denom, Part #Computing securely with coorect degree polynomial: x.tranpose * y def vecVec(x,y): p = len(x) s = secretsharing(0,t,n) for i in range(p): s+= mult(x[i], y[i]) return s #Computing securely with correct degree polynomial: x.transpose * A def vecMat(x,A): p = len(x) res = [] for i in range(p): s = secretsharing(0,t,n) for j in range(p): s+= mult(x[j], A[i+j*p]) res.append(s) return res #Computing securely with correct degree polynomial: A * x def matVec(A,x): res = [] p = len(x) for i in range(p): s = secretsharing(0,t,n) for j in range(p): s += mult(x[j], A[j+(i*p)]) res.append(s) return res #Computing securely with correct degree polynomial: x * x.transpose def VecVec(x,y): res = [] for i in range(len(y)): for j in range(len(x)): res.append(mult(y[i],x[j])) return res 
       
########################################################################################### # This cells computes SECURELY the recursice least square equations using the functions # # from cell 1 and 2 (thus cell 1 and 2 must run before attempting to run this cell). # ########################################################################################### np.random.seed(2) obs = 70 n = 4 # Numbers of parties t = 1 # Corrupted parties p = 2 # Dimension (for one dimensional system, p = 2) C = 2**7 # Scaling constant P = C*np.eye(p) # Initialize P Ps = [] # List for holding shares of P l = 13 # Maximum number of bits for any secret nominator and denominator k=l**2+l # Scaling for 1/d #listObs = [obs*[1]] # 1 dimensional system listObs = [] # Multi dimensional system #params = [2] # 1 dimensional system params = [] # Multi dimensional system for j in range(p): # For 1 dimensional system write p-1. listObs.append(np.random.randint(0,5,obs)) #Observations of input params.append(np.random.randint(1,5,1)) xobs = np.array(listObs).transpose() params = np.array(params).flatten() print 'True parameters: ', params print pretty_print(html('<u> Time | Estimate </u>')) y = [sum(i * params) for i in xobs] #Observations of output noice = np.random.normal(0,2,obs) #Noice y = y+noice #Add noice y = np.array(y,dtype = int) #Truncate Y #The P matrix, will be a list of the secret entries. --> Ps = [P[1,1], P[1,2], P[2,1], P[2,2]]. for i in range(p): for j in range(p): Ps.append(secretsharing(P[i][j],t,n)) w = [secretsharing(0,t,n) for i in range(p)] #Initialization of the parameter estimate wOPEN = [] for i in range(obs): x = xobs[i] #Observation vector for 3 X = [list(secretsharing(j,t,n)) for j in x] #Secret sharing of observations denom, Ppart = estP(F,Ps, X,C,t,n) # TJECK FOR WRAP-AROUND ZERO b = [lessthan(bitDec(-Ppart[j],2*l),bitDec(Ppart[j],2*l)) for j in range(p**2)] #REMOVE SIGN for j in range(p**2): Ppart[j] = mult(b[j], -Ppart[j]) + mult(1-b[j], Ppart[j]) #PERFORM DIVISION P2 = [integerDiv(F, denom, np.array(j), n, t, l, k) for j in Ppart] #ADD SIGN for j in range(p**2): P2[j] = mult(b[j], -P2[j]) + mult(1-b[j], P2[j]) Ps = [h-j for h,j in zip(Ps, P2)] y2 = secretsharing(y[i],t,n) w = ILS(F, Ps, X, w, y2, C,t,n,i) w2 = np.copy(w) w_open = [] #REVEAL (FLOATING POINT) RESULT WITH ADDED SIGN #TJECK WRAP-AROUND ZERO b2 = [less(-j, j) for j in w2] for j in w2: w_open.append(np.array(rec(j),dtype = float)/C) for j in range(p): if b2[j]: w_open[j] = -w_open[j] w_open = np.array(w_open) print ' {:>2} | {} '.format(i, w_open) #print i'th estimate. wOPEN.append(w_open) print '---------------------------------------' e = [] for j in range(obs): e.append(sum((params - wOPEN[j])**2) * 1./p) #print '({},{:.10f})'.format(j, e[j]) #PLOTTING PURPOSE list_plot(e) 
       
True parameters:  [3 2]

 Time |       Estimate           
   0  | [ 0.  2.] 
   1  | [ 0.  2.] 
   2  | [ 2.7421875  2.2109375] 
   3  | [ 2.8984375  2.3671875] 
True parameters:  [3 2]

 Time |       Estimate           
   0  | [ 0.  2.] 
   1  | [ 0.  2.] 
   2  | [ 2.7421875  2.2109375] 
   3  | [ 2.8984375  2.3671875] 
########################################################################################### # This cells computes NON-SECURELY the recursice least square equations. # # Thus FLOATING-point numbers are used in these equations. Can be used as a comparison # # for the results from cell 3. # ########################################################################################### import numpy as np np.random.seed(2) obs = 70 p = 2 #listObs = [obs*[1]] #1 dimensional system listObs = [] #Multi dimensional system #params = [2] #1 dimensional system params = [] #Multi dimensional system for j in range(p): listObs.append(np.random.randint(0,5,obs)) #Observations of input params.append(np.random.randint(1,5,1)) xobs = np.array(listObs).transpose() params = np.array(params).flatten() print 'True parameters: ', params print pretty_print(html('<u> Time | Estimate </u>')) y = [sum(i * params) for i in xobs] #Observations of output noice = np.random.normal(0,2,obs) #Noice y = y+noice #Add noice P = matrix.identity(p) w = matrix(p,1,p*[0]) wOPEN = [] for i in range(obs): xdut = [int(j) for j in xobs[i]] #x = matrix(p,1,[1, x1[i]]) #Observation vector for 1 dimension x = matrix(p,1, xdut) #Observation vector for multi dimensions denom = (1 + x.transpose()* P * x) P = P - matrix(denom[0][0]**(-1) *( P * x * x.transpose() * P)) g = P * x e = y[i] - x.transpose() * w w = matrix(w + g * e) w3 = np.array(w) print ' {:>2} | {} '.format(i, w3.flatten()) #print i'th estimate. wOPEN.append(w3) print '---------------------------------------' params = params.reshape(p,1) e = [] for j in range(obs): e.append(sum((params - wOPEN[j])**2) * 1./p) #print '({},{:.10f})'.format(j, e[j]) list_plot(e) 
       
True parameters:  [3 2]

 Time |       Estimate           
   0  | [ 0.          2.11762403] 
   1  | [ 0.          2.15115503] 
   2  | [ 2.53531671  2.2718844 ] 
   3  | [ 2.69172395  2.4152577 ] 
   4  | [ 2.39316875  2.48128433] 
   5  | [ 2.39316875  2.48128433] 
   6  | [ 2.44977158  2.48128433] 
   7  | [ 2.48587672  2.40907404] 
   8  | [ 2.71556808  2.39515335] 
   9  | [ 2.74861134  2.37725492] 
  10  | [ 3.09387726  2.19023588] 
  11  | [ 3.04924542  2.18632081] 
  12  | [ 3.03041738  2.18466922] 
  13  | [ 3.19065131  2.07250547] 
  14  | [ 3.10178173  2.13471419] 
  15  | [ 3.10855541  2.14750892] 
  16  | [ 3.1004336   2.13216774] 
  17  | [ 3.0185982   2.16764145] 
  18  | [ 3.00469424  2.13594042] 
  19  | [ 3.02555767  2.08607534] 
  20  | [ 3.01797773  2.11308749] 
  21  | [ 2.98912791  2.10627893] 
  22  | [ 2.98590143  2.10134958] 
  23  | [ 2.98644209  2.10055033] 
  24  | [ 2.95219204  2.04963808] 
  25  | [ 2.98043229  2.04718748] 
  26  | [ 2.9763524   2.05483516] 
  27  | [ 3.03582114  2.00883724] 
  28  | [ 3.03582114  2.00883724] 
  29  | [ 3.06328571  1.99894009] 
  30  | [ 3.05771222  1.98869692] 
  31  | [ 3.07531218  1.94712221] 
  32  | [ 3.07198832  1.94693335] 
  33  | [ 3.06926044  1.94251678] 
  34  | [ 3.06942524  1.94278361] 
  35  | [ 3.06882634  1.943598  ] 
  36  | [ 3.02323153  1.96365971] 
  37  | [ 3.01703996  1.95334044] 
  38  | [ 3.016944   1.9531805] 
  39  | [ 2.96738468  2.02043958] 
  40  | [ 2.98300259  2.00804641] 
  41  | [ 2.97150409  2.02371019] 
  42  | [ 3.0126644   1.99118027] 
  43  | [ 3.04172681  1.94303209] 
  44  | [ 3.0459625   1.93526665] 
  45  | [ 3.04501986  1.93432401] 
  46  | [ 3.06667965  1.93459252] 
  47  | [ 3.02997219  1.96260952] 
  48  | [ 3.05994971  1.93972915] 
  49  | [ 3.0466886   1.95833155] 
  50  | [ 3.06731382  1.95638051] 
  51  | [ 3.08440893  1.93235496] 
  52  | [ 3.07557177  1.92215824] 
  53  | [ 3.07342069  1.92292406] 
  54  | [ 3.06578026  1.93364748] 
  55  | [ 3.08295452  1.95201995] 
  56  | [ 3.08025136  1.97154279] 
  57  | [ 3.08413803  1.96625691] 
  58  | [ 3.1180351  1.9413326] 
  59  | [ 3.14417511  1.94270388] 
  60  | [ 3.13931347  1.94635359] 
  61  | [ 3.15150042  1.9602215 ] 
  62  | [ 3.15435226  1.95632785] 
  63  | [ 3.15148839  1.97345079] 
  64  | [ 3.17948488  1.92235569] 
  65  | [ 3.12452563  1.96317422] 
  66  | [ 3.14341512  1.93679641] 
  67  | [ 3.115288    1.98546921] 
  68  | [ 3.11240871  1.99374303] 
  69  | [ 3.11558913  1.99678147] 
---------------------------------------
True parameters:  [3 2]

 Time |       Estimate           
   0  | [ 0.          2.11762403] 
   1  | [ 0.          2.15115503] 
   2  | [ 2.53531671  2.2718844 ] 
   3  | [ 2.69172395  2.4152577 ] 
   4  | [ 2.39316875  2.48128433] 
   5  | [ 2.39316875  2.48128433] 
   6  | [ 2.44977158  2.48128433] 
   7  | [ 2.48587672  2.40907404] 
   8  | [ 2.71556808  2.39515335] 
   9  | [ 2.74861134  2.37725492] 
  10  | [ 3.09387726  2.19023588] 
  11  | [ 3.04924542  2.18632081] 
  12  | [ 3.03041738  2.18466922] 
  13  | [ 3.19065131  2.07250547] 
  14  | [ 3.10178173  2.13471419] 
  15  | [ 3.10855541  2.14750892] 
  16  | [ 3.1004336   2.13216774] 
  17  | [ 3.0185982   2.16764145] 
  18  | [ 3.00469424  2.13594042] 
  19  | [ 3.02555767  2.08607534] 
  20  | [ 3.01797773  2.11308749] 
  21  | [ 2.98912791  2.10627893] 
  22  | [ 2.98590143  2.10134958] 
  23  | [ 2.98644209  2.10055033] 
  24  | [ 2.95219204  2.04963808] 
  25  | [ 2.98043229  2.04718748] 
  26  | [ 2.9763524   2.05483516] 
  27  | [ 3.03582114  2.00883724] 
  28  | [ 3.03582114  2.00883724] 
  29  | [ 3.06328571  1.99894009] 
  30  | [ 3.05771222  1.98869692] 
  31  | [ 3.07531218  1.94712221] 
  32  | [ 3.07198832  1.94693335] 
  33  | [ 3.06926044  1.94251678] 
  34  | [ 3.06942524  1.94278361] 
  35  | [ 3.06882634  1.943598  ] 
  36  | [ 3.02323153  1.96365971] 
  37  | [ 3.01703996  1.95334044] 
  38  | [ 3.016944   1.9531805] 
  39  | [ 2.96738468  2.02043958] 
  40  | [ 2.98300259  2.00804641] 
  41  | [ 2.97150409  2.02371019] 
  42  | [ 3.0126644   1.99118027] 
  43  | [ 3.04172681  1.94303209] 
  44  | [ 3.0459625   1.93526665] 
  45  | [ 3.04501986  1.93432401] 
  46  | [ 3.06667965  1.93459252] 
  47  | [ 3.02997219  1.96260952] 
  48  | [ 3.05994971  1.93972915] 
  49  | [ 3.0466886   1.95833155] 
  50  | [ 3.06731382  1.95638051] 
  51  | [ 3.08440893  1.93235496] 
  52  | [ 3.07557177  1.92215824] 
  53  | [ 3.07342069  1.92292406] 
  54  | [ 3.06578026  1.93364748] 
  55  | [ 3.08295452  1.95201995] 
  56  | [ 3.08025136  1.97154279] 
  57  | [ 3.08413803  1.96625691] 
  58  | [ 3.1180351  1.9413326] 
  59  | [ 3.14417511  1.94270388] 
  60  | [ 3.13931347  1.94635359] 
  61  | [ 3.15150042  1.9602215 ] 
  62  | [ 3.15435226  1.95632785] 
  63  | [ 3.15148839  1.97345079] 
  64  | [ 3.17948488  1.92235569] 
  65  | [ 3.12452563  1.96317422] 
  66  | [ 3.14341512  1.93679641] 
  67  | [ 3.115288    1.98546921] 
  68  | [ 3.11240871  1.99374303] 
  69  | [ 3.11558913  1.99678147] 
---------------------------------------