Privacy Preserving Recursive Least Squares

17 days ago by ktjell12@student.aau.dk

##################################################################################################### # This cell implementes most of the SMPC funtions, used in the paper: # # Privacy Preserving Recursive Least Squares Solutions by Katrine T., Rafael W., Ignacio C. # # # # To name a few: # # - Shamir's Secret sharing scheme # # - Multiplication protocol # # - Integer division # # - Comparison protocol # ##################################################################################################### import numpy as np import matplotlib.pylab as plt import csv import sys #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 = x for j in range(1,t+1): s = F(s) + F(c[j-1]) * F(i)**F(j) shares.append(s) return np.array(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 += F(x[i] * y[i]) return res #Create a Beaver's triplet. def triplet(n, t): r = basispoly(n) matrixA = np.zeros((n,n)) matrixB = np.zeros((n,n)) 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 = np.array(a)*np.array(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+= np.array(r[i]) * np.array(hs[i]) return np.array(s) #Integer division with public constant #def divPub(x, k, l): # Should return a random bit unknown to all parties. Does not seem to work for huge field. 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 print rec(r_shares) rsq_shares = mult(r_shares,r_shares) print rec(rsq_shares) #salpha, sbeta, sgamma = triplet(n, t) #d = dot(F, basis, [F(i-j) for (i,j) in zip(r_shares, salpha)]) #e = dot(F, basis, [F(i-j) for (i,j) in zip(r_shares, sbeta)]) #rsq_shares = [F(d*e + d*i + e*j + k) for (i,j,k) in zip(sbeta, salpha, sgamma)] rsq = rec(rsq_shares) #r**2 is revealed r_prime = sqrt(rsq) print 'rp', r_prime #r = np.random.choice(2,1)[0] return (F(1/2) * (r_prime**(-1) * np.array(r_shares) + 1)) #Return both sharings of a random number and the bitwise sharing of the same random number. def randomBitNumber(F, n, t, l): sharings = [] for i in range(l): sharings.append(randomBit(F,n,t)) sharings_reverse = sharings[::-1] r_shares = [] for i in range(n): s = 0 for j in range(l): s += 2**j * rec(sharings_reverse[j]) r_shares.append(s) return sharings, r_shares #Computing the inverse of a field element. def inverse(x, r): w_s = [i*j for i,j in zip(x,r)] w = rec(w_s) return w**(-1) * np.array(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 = np.array(a1) - np.array(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+= a[i] * 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)): f.append(OR(e[0:i+1])) g = [f[0]] for i in range(0,len(f)-1): g.append(f[i+1] - f[i]) h = [] for i in range(len(g)): h.append(mult(g[i], b[i])) y = [] for i in range(len(h[0])): s = 0 for j in range(len(h)): s += h[j][i] y.append(s) return y #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 #Computing PRE-OR more efficiently. def pre(F, d, r, rb, n, t, l): ### ADD RANDOMNESS basis = basispoly(n) rb.insert(0, list(secretsharing(0,t,n))) es = np.array(d) + np.array(r) e = rec(es) ebit = list(bin(e)[2:].rjust(l+1,'0')) ebs = [] for b in ebit: ebs.append(secretsharing(b,t,n)) ### Create H ep = [xor(i,j) for i,j in zip(ebs,rb)] E = [] for i in range(len(ep)): E.append(OR(ep[0:i+1])) h = [i - mult(i,j) for i,j in zip (rb, ebs)] hp = [] for i in range(len(h)): if i ==0: Es = 0 else: Es = E[i-1] hp.append(h[i] + 1 - Es) H = [] for i in range(len(hp)): H.append(AND(hp[0:i+1])) ### 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.append(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_] D = D[::-1] 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 y #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(F, n, t, l) rr = rec(r) D = pre(F, d, r, rbits, n, t, l) shift = [1+i for i in D] shiftInv = inverse(shift, r) a1 = [i-j for i,j in zip(shift, d)] p = mult(a1,shiftInv) sum_s = [] prod = secretsharing(1,t,n) s = np.array(n*[0]) for i in range(1,l): prod = mult(prod, p) s = s + np.array(prod) s = [1+i for i in s] a_hat = [2**k*i for i in shiftInv] a_hat = mult(a_hat,s) q_hat = mult(nom,a_hat) q = down(F, q_hat, 2**k,t,n) r = nom - mult(d,q) dd = rec(d+d) rd = rec(r+d) ep = int(dd <= rd) em = int(rec(d) > rd) 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 = 13 # bit length of some secrets r = basispoly(n) # Lagrange basispolynomial coefficients print rec(randomBit(n,t)) 
       
625198650575593907948139785385919338907056866949899424242338243422238322\
81
634581691528657710873525423649129981338764148962396127626508371755731539\
231
rp
sqrt63458169152865771087352542364912998133876414896239612762650837175573\
1539231
Traceback (click to the left of this block for traceback)
...
ValueError: element is not in the prime field
62519865057559390794813978538591933890705686694989942424233824342223832281
634581691528657710873525423649129981338764148962396127626508371755731539231
rp sqrt634581691528657710873525423649129981338764148962396127626508371755731539231
Traceback (most recent call last):    #    - Multiplication protocol                                                                      #
  File "", line 1, in <module>
    
  File "/tmp/tmpxi0y4I/___code___.py", line 316, in <module>
    exec compile(u'print rec(randomBit(n,t))
  File "", line 1, in <module>
    
  File "/tmp/tmpxi0y4I/___code___.py", line 46, in rec
    res += F(x[i] * y[i])
  File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873)
  File "sage/rings/finite_rings/hom_prime_finite_field.pyx", line 44, in sage.rings.finite_rings.hom_prime_finite_field.SectionFiniteFieldHomomorphism_prime._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/hom_prime_finite_field.c:3132)
  File "/usr/local/sage-6.9/local/lib/python2.7/site-packages/sage/rings/finite_rings/integer_mod_ring.py", line 1167, in _element_constructor_
    return integer_mod.IntegerMod(self, x)
  File "sage/rings/finite_rings/integer_mod.pyx", line 182, in sage.rings.finite_rings.integer_mod.IntegerMod (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:4160)
  File "sage/rings/finite_rings/integer_mod.pyx", line 1716, in sage.rings.finite_rings.integer_mod.IntegerMod_gmp.__init__ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:20448)
  File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873)
  File "sage/structure/coerce_maps.pyx", line 237, in sage.structure.coerce_maps.NamedConvertMap._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/coerce_maps.c:5944)
  File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 937, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt._integer_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9355)
  File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 920, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt.lift (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9228)
ValueError: element is not in the prime field
random_prime(2**250, None, 2**247) 
       
696205902887547871652159517929547663282608944037200206897421778089495766\
033
696205902887547871652159517929547663282608944037200206897421778089495766033
import numpy as np import matplotlib.pylab as plt import csv import sys #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 = x for j in range(1,t+1): s = F(s) + F(c[j-1]) * F(i)**F(j) shares.append(s) return np.array(shares) g #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 += F(x[i] * y[i]) return res #Create a Beaver's triplet. def triplet(n, t): r = basispoly(n) matrixA = np.zeros((n,n)) matrixB = np.zeros((n,n)) 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 = np.array(a)*np.array(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+= np.array(r[i]) * np.array(hs[i]) return np.array(s) #Integer division with public constant #def divPub(x, k, l): 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 = 13 # bit length of some secrets r = basispoly(n) # Lagrange basispolynomial coefficients 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) * (r_prime**(-1) * np.array(r_shares) + 1)) print rec(randomBit(n,t)) 
       
Traceback (click to the left of this block for traceback)
...
ValueError: element is not in the prime field
Traceback (most recent call last):        shares = []
  File "", line 1, in <module>
    
  File "/tmp/tmpbir8B8/___code___.py", line 97, in <module>
    exec compile(u'print rec(randomBit(n,t))
  File "", line 1, in <module>
    
  File "/tmp/tmpbir8B8/___code___.py", line 36, in rec
    res += F(x[i] * y[i])
  File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873)
  File "sage/rings/finite_rings/hom_prime_finite_field.pyx", line 44, in sage.rings.finite_rings.hom_prime_finite_field.SectionFiniteFieldHomomorphism_prime._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/hom_prime_finite_field.c:3132)
  File "/usr/local/sage-6.9/local/lib/python2.7/site-packages/sage/rings/finite_rings/integer_mod_ring.py", line 1167, in _element_constructor_
    return integer_mod.IntegerMod(self, x)
  File "sage/rings/finite_rings/integer_mod.pyx", line 182, in sage.rings.finite_rings.integer_mod.IntegerMod (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:4160)
  File "sage/rings/finite_rings/integer_mod.pyx", line 1716, in sage.rings.finite_rings.integer_mod.IntegerMod_gmp.__init__ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:20448)
  File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873)
  File "sage/structure/coerce_maps.pyx", line 237, in sage.structure.coerce_maps.NamedConvertMap._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/coerce_maps.c:5944)
  File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 937, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt._integer_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9355)
  File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 920, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt.lift (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9228)
ValueError: element is not in the prime field
 
       
b########################################################################### #This cells defines functions that are more specifically used to # #compute the recursive least squares equations securely. # ########################################################################### #Should and could, be implemented securely. def less(a,b): return int( rec(a) < rec(b) ) 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) 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 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 (COULD and SHOULD be implemented securely) b = [less(Ppart[j]*-1,Ppart[j]) for j in range(p**2)] #REMOVE SIGN for j in range(p**2): if b[j]: Ppart[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): if b[j]: P2[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 = [] #TJECK WRAP-AROUND ZERO b2 = [less(j*-1, j) for j in w2] for j in range(p): if b2[j]: w2[j] = -w2[j] #REVEAL (FLOATING POINT) RESULT WITH ADDED SIGN 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] 
   4  | [ 2.5234375  2.4609375] 
   5  | [ 2.5234375  2.4609375] 
   6  | [ 2.5234375  2.4609375] 
   7  | [ 2.5703125  2.3671875] 
   8  | [ 2.921875   2.3671875] 
   9  | [ 2.921875   2.3671875] 
  10  | [ 3.390625   2.0859375] 
  11  | [ 3.3046875  2.0625   ] 
  12  | [ 3.21875    2.0390625] 
  13  | [ 3.40625    1.8984375] 
  14  | [ 3.21875    2.0234375] 
  15  | [ 3.234375   2.0703125] 
  16  | [ 3.234375   2.0703125] 
  17  | [ 3.0703125  2.09375  ] 
  18  | [ 3.0390625  2.       ] 
  19  | [ 3.0625    1.921875] 
  20  | [ 3.03125   2.109375] 
  21  | [ 2.984375  2.078125] 
  22  | [ 2.9609375  2.03125  ] 
  23  | [ 2.9609375  2.03125  ] 
  24  | [ 2.8671875  1.84375  ] 
  25  | [ 2.9453125  1.875    ] 
  26  | [ 2.921875  1.953125] 
  27  | [ 3.0625    1.859375] 
  28  | [ 3.0625    1.859375] 
  29  | [ 3.125     1.859375] 
  30  | [ 3.125     1.859375] 
  31  | [ 3.140625  1.765625] 
  32  | [ 3.140625  1.765625] 
  33  | [ 3.140625  1.765625] 
  34  | [ 3.1484375  1.78125  ] 
  35  | [ 3.1328125  1.8125   ] 
  36  | [ 3.0078125  1.8125   ] 
  37  | [ 2.9921875  1.78125  ] 
  38  | [ 2.9921875  1.78125  ] 
  39  | [ 2.8515625  2.0625   ] 
  40  | [ 2.8984375  2.03125  ] 
  41  | [ 2.8828125  2.0625   ] 
  42  | [ 3.0234375  1.96875  ] 
  43  | [ 3.0625    1.890625] 
  44  | [ 3.0625    1.890625] 
  45  | [ 3.0625    1.890625] 
  46  | [ 3.15625  1.90625] 
  47  | [ 3.015625  2.      ] 
  48  | [ 3.0859375  1.953125 ] 
  49  | [ 3.0390625  2.0234375] 
  50  | [ 3.1171875  2.0234375] 
  51  | [ 3.1796875  1.9296875] 
  52  | [ 3.1328125  1.8828125] 
  53  | [ 3.1015625  1.890625 ] 
  54  | [ 3.0546875  1.9609375] 
  55  | [ 3.1171875  2.0234375] 
  56  | [ 3.1171875  2.1015625] 
  57  | [ 3.1328125  2.078125 ] 
  58  | [ 3.2265625  2.015625 ] 
  59  | [ 3.3203125  2.03125  ] 
  60  | [ 3.296875  2.046875] 
  61  | [ 3.34375  2.09375] 
  62  | [ 3.359375   2.0703125] 
  63  | [ 3.359375   2.1484375] 
  64  | [ 3.453125   1.9296875] 
  65  | [ 3.171875   2.1171875] 
  66  | [ 3.234375   2.0234375] 
  67  | [ 3.15625    2.1796875] 
  68  | [ 3.1484375  2.2109375] 
  69  | [ 3.1640625  2.2265625] 
---------------------------------------
True parameters:  [3 2]

 Time |       Estimate           
   0  | [ 0.  2.] 
   1  | [ 0.  2.] 
   2  | [ 2.7421875  2.2109375] 
   3  | [ 2.8984375  2.3671875] 
   4  | [ 2.5234375  2.4609375] 
   5  | [ 2.5234375  2.4609375] 
   6  | [ 2.5234375  2.4609375] 
   7  | [ 2.5703125  2.3671875] 
   8  | [ 2.921875   2.3671875] 
   9  | [ 2.921875   2.3671875] 
  10  | [ 3.390625   2.0859375] 
  11  | [ 3.3046875  2.0625   ] 
  12  | [ 3.21875    2.0390625] 
  13  | [ 3.40625    1.8984375] 
  14  | [ 3.21875    2.0234375] 
  15  | [ 3.234375   2.0703125] 
  16  | [ 3.234375   2.0703125] 
  17  | [ 3.0703125  2.09375  ] 
  18  | [ 3.0390625  2.       ] 
  19  | [ 3.0625    1.921875] 
  20  | [ 3.03125   2.109375] 
  21  | [ 2.984375  2.078125] 
  22  | [ 2.9609375  2.03125  ] 
  23  | [ 2.9609375  2.03125  ] 
  24  | [ 2.8671875  1.84375  ] 
  25  | [ 2.9453125  1.875    ] 
  26  | [ 2.921875  1.953125] 
  27  | [ 3.0625    1.859375] 
  28  | [ 3.0625    1.859375] 
  29  | [ 3.125     1.859375] 
  30  | [ 3.125     1.859375] 
  31  | [ 3.140625  1.765625] 
  32  | [ 3.140625  1.765625] 
  33  | [ 3.140625  1.765625] 
  34  | [ 3.1484375  1.78125  ] 
  35  | [ 3.1328125  1.8125   ] 
  36  | [ 3.0078125  1.8125   ] 
  37  | [ 2.9921875  1.78125  ] 
  38  | [ 2.9921875  1.78125  ] 
  39  | [ 2.8515625  2.0625   ] 
  40  | [ 2.8984375  2.03125  ] 
  41  | [ 2.8828125  2.0625   ] 
  42  | [ 3.0234375  1.96875  ] 
  43  | [ 3.0625    1.890625] 
  44  | [ 3.0625    1.890625] 
  45  | [ 3.0625    1.890625] 
  46  | [ 3.15625  1.90625] 
  47  | [ 3.015625  2.      ] 
  48  | [ 3.0859375  1.953125 ] 
  49  | [ 3.0390625  2.0234375] 
  50  | [ 3.1171875  2.0234375] 
  51  | [ 3.1796875  1.9296875] 
  52  | [ 3.1328125  1.8828125] 
  53  | [ 3.1015625  1.890625 ] 
  54  | [ 3.0546875  1.9609375] 
  55  | [ 3.1171875  2.0234375] 
  56  | [ 3.1171875  2.1015625] 
  57  | [ 3.1328125  2.078125 ] 
  58  | [ 3.2265625  2.015625 ] 
  59  | [ 3.3203125  2.03125  ] 
  60  | [ 3.296875  2.046875] 
  61  | [ 3.34375  2.09375] 
  62  | [ 3.359375   2.0703125] 
  63  | [ 3.359375   2.1484375] 
  64  | [ 3.453125   1.9296875] 
  65  | [ 3.171875   2.1171875] 
  66  | [ 3.234375   2.0234375] 
  67  | [ 3.15625    2.1796875] 
  68  | [ 3.1484375  2.2109375] 
  69  | [ 3.1640625  2.2265625] 
---------------------------------------
########################################################################################### # 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, w_open) #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  | [ 1.875      2.3671875] 
   1  | [ 1.875      2.3671875] 
   2  | [ 1.875      2.3671875] 
   3  | [ 1.875      2.3671875] 
   4  | [ 1.875      2.3671875] 
   5  | [ 1.875      2.3671875] 
   6  | [ 1.875      2.3671875] 
   7  | [ 1.875      2.3671875] 
   8  | [ 1.875      2.3671875] 
   9  | [ 1.875      2.3671875] 
  10  | [ 1.875      2.3671875] 
  11  | [ 1.875      2.3671875] 
  12  | [ 1.875      2.3671875] 
  13  | [ 1.875      2.3671875] 
  14  | [ 1.875      2.3671875] 
  15  | [ 1.875      2.3671875] 
  16  | [ 1.875      2.3671875] 
  17  | [ 1.875      2.3671875] 
  18  | [ 1.875      2.3671875] 
  19  | [ 1.875      2.3671875] 
  20  | [ 1.875      2.3671875] 
  21  | [ 1.875      2.3671875] 
  22  | [ 1.875      2.3671875] 
  23  | [ 1.875      2.3671875] 
  24  | [ 1.875      2.3671875] 
  25  | [ 1.875      2.3671875] 
  26  | [ 1.875      2.3671875] 
  27  | [ 1.875      2.3671875] 
  28  | [ 1.875      2.3671875] 
  29  | [ 1.875      2.3671875] 
  30  | [ 1.875      2.3671875] 
  31  | [ 1.875      2.3671875] 
  32  | [ 1.875      2.3671875] 
  33  | [ 1.875      2.3671875] 
  34  | [ 1.875      2.3671875] 
  35  | [ 1.875      2.3671875] 
  36  | [ 1.875      2.3671875] 
  37  | [ 1.875      2.3671875] 
  38  | [ 1.875      2.3671875] 
  39  | [ 1.875      2.3671875] 
  40  | [ 1.875      2.3671875] 
  41  | [ 1.875      2.3671875] 
  42  | [ 1.875      2.3671875] 
  43  | [ 1.875      2.3671875] 
  44  | [ 1.875      2.3671875] 
  45  | [ 1.875      2.3671875] 
  46  | [ 1.875      2.3671875] 
  47  | [ 1.875      2.3671875] 
  48  | [ 1.875      2.3671875] 
  49  | [ 1.875      2.3671875] 
  50  | [ 1.875      2.3671875] 
  51  | [ 1.875      2.3671875] 
  52  | [ 1.875      2.3671875] 
  53  | [ 1.875      2.3671875] 
  54  | [ 1.875      2.3671875] 
  55  | [ 1.875      2.3671875] 
  56  | [ 1.875      2.3671875] 
  57  | [ 1.875      2.3671875] 
  58  | [ 1.875      2.3671875] 
  59  | [ 1.875      2.3671875] 
  60  | [ 1.875      2.3671875] 
  61  | [ 1.875      2.3671875] 
  62  | [ 1.875      2.3671875] 
  63  | [ 1.875      2.3671875] 
  64  | [ 1.875      2.3671875] 
  65  | [ 1.875      2.3671875] 
  66  | [ 1.875      2.3671875] 
  67  | [ 1.875      2.3671875] 
  68  | [ 1.875      2.3671875] 
  69  | [ 1.875      2.3671875] 
---------------------------------------
True parameters:  [3 2]

 Time |       Estimate           
   0  | [ 1.875      2.3671875] 
   1  | [ 1.875      2.3671875] 
   2  | [ 1.875      2.3671875] 
   3  | [ 1.875      2.3671875] 
   4  | [ 1.875      2.3671875] 
   5  | [ 1.875      2.3671875] 
   6  | [ 1.875      2.3671875] 
   7  | [ 1.875      2.3671875] 
   8  | [ 1.875      2.3671875] 
   9  | [ 1.875      2.3671875] 
  10  | [ 1.875      2.3671875] 
  11  | [ 1.875      2.3671875] 
  12  | [ 1.875      2.3671875] 
  13  | [ 1.875      2.3671875] 
  14  | [ 1.875      2.3671875] 
  15  | [ 1.875      2.3671875] 
  16  | [ 1.875      2.3671875] 
  17  | [ 1.875      2.3671875] 
  18  | [ 1.875      2.3671875] 
  19  | [ 1.875      2.3671875] 
  20  | [ 1.875      2.3671875] 
  21  | [ 1.875      2.3671875] 
  22  | [ 1.875      2.3671875] 
  23  | [ 1.875      2.3671875] 
  24  | [ 1.875      2.3671875] 
  25  | [ 1.875      2.3671875] 
  26  | [ 1.875      2.3671875] 
  27  | [ 1.875      2.3671875] 
  28  | [ 1.875      2.3671875] 
  29  | [ 1.875      2.3671875] 
  30  | [ 1.875      2.3671875] 
  31  | [ 1.875      2.3671875] 
  32  | [ 1.875      2.3671875] 
  33  | [ 1.875      2.3671875] 
  34  | [ 1.875      2.3671875] 
  35  | [ 1.875      2.3671875] 
  36  | [ 1.875      2.3671875] 
  37  | [ 1.875      2.3671875] 
  38  | [ 1.875      2.3671875] 
  39  | [ 1.875      2.3671875] 
  40  | [ 1.875      2.3671875] 
  41  | [ 1.875      2.3671875] 
  42  | [ 1.875      2.3671875] 
  43  | [ 1.875      2.3671875] 
  44  | [ 1.875      2.3671875] 
  45  | [ 1.875      2.3671875] 
  46  | [ 1.875      2.3671875] 
  47  | [ 1.875      2.3671875] 
  48  | [ 1.875      2.3671875] 
  49  | [ 1.875      2.3671875] 
  50  | [ 1.875      2.3671875] 
  51  | [ 1.875      2.3671875] 
  52  | [ 1.875      2.3671875] 
  53  | [ 1.875      2.3671875] 
  54  | [ 1.875      2.3671875] 
  55  | [ 1.875      2.3671875] 
  56  | [ 1.875      2.3671875] 
  57  | [ 1.875      2.3671875] 
  58  | [ 1.875      2.3671875] 
  59  | [ 1.875      2.3671875] 
  60  | [ 1.875      2.3671875] 
  61  | [ 1.875      2.3671875] 
  62  | [ 1.875      2.3671875] 
  63  | [ 1.875      2.3671875] 
  64  | [ 1.875      2.3671875] 
  65  | [ 1.875      2.3671875] 
  66  | [ 1.875      2.3671875] 
  67  | [ 1.875      2.3671875] 
  68  | [ 1.875      2.3671875] 
  69  | [ 1.875      2.3671875] 
---------------------------------------