1import numpy 2from algopy import CGraph, Function, UTPM, dot, qr, qr_full, eigh, inv, solve, zeros 3 4def eval_covariance_matrix_naive(J1, J2): 5 M,N = J1.shape 6 K,N = J2.shape 7 tmp = zeros((N+K, N+K), dtype=J1) 8 tmp[:N,:N] = dot(J1.T,J1) 9 tmp[:N,N:] = J2.T 10 tmp[N:,:N] = J2 11 return inv(tmp)[:N,:N] 12 13def eval_covariance_matrix_qr(J1, J2): 14 M,N = J1.shape 15 K,N = J2.shape 16 Q,R = qr_full(J2.T) 17 Q2 = Q[:,K:].T 18 J1_tilde = dot(J1,Q2.T) 19 Q,R = qr(J1_tilde) 20 V = solve(R.T, Q2) 21 return dot(V.T,V) 22 23 24# dimensions of the involved matrices 25D,P,M,N,K,Nx = 2,1,5,3,1,1 26 27# trace the function evaluation of METHOD 1: nullspace method 28cg1 = CGraph() 29J1 = Function(UTPM(numpy.random.rand(*(D,P,M,N)))) 30J2 = Function(UTPM(numpy.random.rand(*(D,P,K,N)))) 31C = eval_covariance_matrix_qr(J1, J2) 32y = C[0,0] 33cg1.trace_off() 34cg1.independentFunctionList = [J1, J2] 35cg1.dependentFunctionList = [y] 36print('covariance matrix: C =\n',C) 37 38# trace the function evaluation of METHOD 2: naive method (potentially numerically unstable) 39cg2 = CGraph() 40J1 = Function(J1.x) 41J2 = Function(J2.x) 42C2 = eval_covariance_matrix_naive(J1, J2) 43y = C2[0,0] 44cg2.trace_off() 45cg2.independentFunctionList = [J1, J2] 46cg2.dependentFunctionList = [y] 47print('covariance matrix: C =\n',C2) 48 49# check that both algorithms returns the same result 50print('difference between naive and nullspace method:\n',C - C2) 51 52# compute the gradient for another value of J1 and J2 53J1 = numpy.random.rand(*(M,N)) 54J2 = numpy.random.rand(*(K,N)) 55 56g1 = cg1.gradient([J1,J2]) 57g2 = cg2.gradient([J1,J2]) 58 59print('naive approach: dy/dJ1 = ', g1[0]) 60print('naive approach: dy/dJ2 = ', g1[1]) 61 62print('nullspace approach: dy/dJ1 = ', g2[0]) 63print('nullspace approach: dy/dJ2 = ', g2[1]) 64 65 66 67 68 69