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