def exKxz_pairwise(self, Z, Xmu, Xcov):
"""
<x_t K_{x_{t-1}, Z}>_q_{x_{t-1:t}}
:param Z: MxD inducing inputs
:param Xmu: X mean (N+1xD)
:param Xcov: 2x(N+1)xDxD
:return: NxMxD
"""
msg_input_shape = "Currently cannot handle slicing in exKxz_pairwise."
assert_input_shape = tf.assert_equal(tf.shape(Xmu)[1], self.input_dim, message=msg_input_shape)
assert_cov_shape = tf.assert_equal(tf.shape(Xmu), tf.shape(Xcov)[1:3], name="assert_Xmu_Xcov_shape")
with tf.control_dependencies([assert_input_shape, assert_cov_shape]):
Xmu = tf.identity(Xmu)
N = tf.shape(Xmu)[0] - 1
D = tf.shape(Xmu)[1]
Xsigmb = tf.slice(Xcov, [0, 0, 0, 0], tf.stack([-1, N, -1, -1]))
Xsigm = Xsigmb[0, :, :, :] # NxDxD
Xsigmc = Xsigmb[1, :, :, :] # NxDxD
Xmum = tf.slice(Xmu, [0, 0], tf.stack([N, -1]))
Xmup = Xmu[1:, :]
lengthscales = self.lengthscales if self.ARD else tf.zeros((D,), dtype=settings.float_type) + self.lengthscales
scalemat = tf.expand_dims(tf.matrix_diag(lengthscales ** 2.0), 0) + Xsigm # NxDxD
det = tf.matrix_determinant(
tf.expand_dims(tf.eye(tf.shape(Xmu)[1], dtype=settings.float_type), 0) +
tf.reshape(lengthscales ** -2.0, (1, 1, -1)) * Xsigm) # N
vec = tf.expand_dims(tf.transpose(Z), 0) - tf.expand_dims(Xmum, 2) # NxDxM
smIvec = tf.matrix_solve(scalemat, vec) # NxDxM
q = tf.reduce_sum(smIvec * vec, [1]) # NxM
addvec = tf.matmul(smIvec, Xsigmc, transpose_a=True) + tf.expand_dims(Xmup, 1) # NxMxD
return self.variance * addvec * tf.reshape(det ** -0.5, (N, 1, 1)) * tf.expand_dims(tf.exp(-0.5 * q), 2)
python类matrix_solve()的实例源码
def exKxz(self, Z, Xmu, Xcov):
"""
It computes the expectation:
<x_t K_{x_t, Z}>_q_{x_t}
:param Z: MxD inducing inputs
:param Xmu: X mean (NxD)
:param Xcov: NxDxD
:return: NxMxD
"""
msg_input_shape = "Currently cannot handle slicing in exKxz."
assert_input_shape = tf.assert_equal(tf.shape(Xmu)[1], self.input_dim, message=msg_input_shape)
assert_cov_shape = tf.assert_equal(tf.shape(Xmu), tf.shape(Xcov)[:2], name="assert_Xmu_Xcov_shape")
with tf.control_dependencies([assert_input_shape, assert_cov_shape]):
Xmu = tf.identity(Xmu)
N = tf.shape(Xmu)[0]
D = tf.shape(Xmu)[1]
lengthscales = self.lengthscales if self.ARD else tf.zeros((D,), dtype=settings.float_type) + self.lengthscales
scalemat = tf.expand_dims(tf.matrix_diag(lengthscales ** 2.0), 0) + Xcov # NxDxD
det = tf.matrix_determinant(
tf.expand_dims(tf.eye(tf.shape(Xmu)[1], dtype=settings.float_type), 0) +
tf.reshape(lengthscales ** -2.0, (1, 1, -1)) * Xcov) # N
vec = tf.expand_dims(tf.transpose(Z), 0) - tf.expand_dims(Xmu, 2) # NxDxM
smIvec = tf.matrix_solve(scalemat, vec) # NxDxM
q = tf.reduce_sum(smIvec * vec, [1]) # NxM
addvec = tf.matmul(smIvec, Xcov, transpose_a=True) + tf.expand_dims(Xmu, 1) # NxMxD
return self.variance * addvec * tf.reshape(det ** -0.5, (N, 1, 1)) * tf.expand_dims(tf.exp(-0.5 * q), 2)
def Linear_RBF_eKxzKzx(self, Ka, Kb, Z, Xmu, Xcov):
Xcov = self._slice_cov(Xcov)
Z, Xmu = self._slice(Z, Xmu)
lin, rbf = (Ka, Kb) if isinstance(Ka, Linear) else (Kb, Ka)
if not isinstance(lin, Linear):
TypeError("{in_lin} is not {linear}".format(in_lin=str(type(lin)), linear=str(Linear)))
if not isinstance(rbf, RBF):
TypeError("{in_rbf} is not {rbf}".format(in_rbf=str(type(rbf)), rbf=str(RBF)))
if lin.ARD or type(lin.active_dims) is not slice or type(rbf.active_dims) is not slice:
raise NotImplementedError("Active dims and/or Linear ARD not implemented. "
"Switching to quadrature.")
D = tf.shape(Xmu)[1]
M = tf.shape(Z)[0]
N = tf.shape(Xmu)[0]
if rbf.ARD:
lengthscales = rbf.lengthscales
else:
lengthscales = tf.zeros((D, ), dtype=settings.float_type) + rbf.lengthscales
lengthscales2 = lengthscales ** 2.0
const = rbf.variance * lin.variance * tf.reduce_prod(lengthscales)
gaussmat = Xcov + tf.matrix_diag(lengthscales2)[None, :, :] # NxDxD
det = tf.matrix_determinant(gaussmat) ** -0.5 # N
cgm = tf.cholesky(gaussmat) # NxDxD
tcgm = tf.tile(cgm[:, None, :, :], [1, M, 1, 1])
vecmin = Z[None, :, :] - Xmu[:, None, :] # NxMxD
d = tf.matrix_triangular_solve(tcgm, vecmin[:, :, :, None]) # NxMxDx1
exp = tf.exp(-0.5 * tf.reduce_sum(d ** 2.0, [2, 3])) # NxM
# exp = tf.Print(exp, [tf.shape(exp)])
vecplus = (Z[None, :, :, None] / lengthscales2[None, None, :, None] +
tf.matrix_solve(Xcov, Xmu[:, :, None])[:, None, :, :]) # NxMxDx1
mean = tf.cholesky_solve(
tcgm, tf.matmul(tf.tile(Xcov[:, None, :, :], [1, M, 1, 1]), vecplus))
mean = mean[:, :, :, 0] * lengthscales2[None, None, :] # NxMxD
a = tf.matmul(tf.tile(Z[None, :, :], [N, 1, 1]),
mean * exp[:, :, None] * det[:, None, None] * const, transpose_b=True)
return a + tf.transpose(a, [0, 2, 1])
def test_solve(self):
with self.test_session():
for batch_shape in [(), (2, 3,)]:
for k in [1, 4]:
operator, mat = self._build_operator_and_mat(batch_shape, k)
# Work with 5 simultaneous systems. 5 is arbitrary.
x = self._rng.randn(*(batch_shape + (k, 5)))
self._compare_results(
expected=tf.matrix_solve(mat, x).eval(), actual=operator.solve(x))
def test_sqrt_solve(self):
# Square roots are not unique, but we should still have
# S^{-T} S^{-1} x = A^{-1} x.
# In our case, we should have S = S^T, so then S^{-1} S^{-1} x = A^{-1} x.
with self.test_session():
for batch_shape in [(), (2, 3,)]:
for k in [1, 4]:
operator, mat = self._build_operator_and_mat(batch_shape, k)
# Work with 5 simultaneous systems. 5 is arbitrary.
x = self._rng.randn(*(batch_shape + (k, 5)))
self._compare_results(
expected=tf.matrix_solve(mat, x).eval(),
actual=operator.sqrt_solve(operator.sqrt_solve(x)))
def testSolve(self):
with self.test_session():
for batch_shape in [(), (2, 3,)]:
for k in [1, 4]:
operator, mat = self._build_operator_and_mat(batch_shape, k)
# Work with 5 simultaneous systems. 5 is arbitrary.
x = self._rng.randn(*(batch_shape + (k, 5)))
self._compare_results(
expected=tf.matrix_solve(mat, x).eval(), actual=operator.solve(x))
def testSqrtSolve(self):
# Square roots are not unique, but we should still have
# S^{-T} S^{-1} x = A^{-1} x.
# In our case, we should have S = S^T, so then S^{-1} S^{-1} x = A^{-1} x.
with self.test_session():
for batch_shape in [(), (2, 3,)]:
for k in [1, 4]:
operator, mat = self._build_operator_and_mat(batch_shape, k)
# Work with 5 simultaneous systems. 5 is arbitrary.
x = self._rng.randn(*(batch_shape + (k, 5)))
self._compare_results(
expected=tf.matrix_solve(mat, x).eval(),
actual=operator.sqrt_solve(operator.sqrt_solve(x)))
def block_CG(A_,B_):
"""
block version of CG. Get solution to matrix equation AX = B, ie
X = A^-1 * B. Will be much faster than Cholesky for large-scale problems.
"""
n = tf.shape(B_)[0]
m = tf.shape(B_)[1]
X = tf.zeros((n,m))
V_ = tf.zeros((n,m))
R = B_
R_ = tf.matrix_set_diag(tf.zeros((n,m)),tf.ones([m]))
#somewhat arbitrary again, may want to check sensitivity
CG_EPS = tf.cast(n/1000,"float")
MAX_ITER = tf.div(n,250) + 3
def cond(i,X,R_,R,V_):
return tf.logical_and(i < MAX_ITER, tf.norm(R) > CG_EPS)
def body(i,X,R_,R,V_):
S = tf.matrix_solve(tf.matmul(tf.transpose(R_),R_),
tf.matmul(tf.transpose(R),R))
V = R + tf.matmul(V_,S)
T = tf.matrix_solve(tf.matmul(tf.transpose(V),tf.matmul(A_,V)),
tf.matmul(tf.transpose(R),R))
X = X + tf.matmul(V,T)
V_ = V
R_ = R
R = R - tf.matmul(A_,tf.matmul(V,T))
return i+1,X,R_,R,V_
i = tf.constant(0)
i,X,_,_,_ = tf.while_loop(cond,body,[i,X,R_,R,V_])
return X
def test_MatrixSolve(self):
t = tf.matrix_solve(*self.random((2, 3, 3, 3), (2, 3, 3, 1)), adjoint=False)
self.check(t)
t = tf.matrix_solve(*self.random((2, 3, 3, 3), (2, 3, 3, 1)), adjoint=True)
self.check(t)
def genPerturbations(opt):
with tf.name_scope("genPerturbations"):
X = np.tile(opt.canon4pts[:,0],[opt.batchSize,1])
Y = np.tile(opt.canon4pts[:,1],[opt.batchSize,1])
dX = tf.random_normal([opt.batchSize,4])*opt.pertScale \
+tf.random_normal([opt.batchSize,1])*opt.transScale
dY = tf.random_normal([opt.batchSize,4])*opt.pertScale \
+tf.random_normal([opt.batchSize,1])*opt.transScale
O = np.zeros([opt.batchSize,4],dtype=np.float32)
I = np.ones([opt.batchSize,4],dtype=np.float32)
# fit warp parameters to generated displacements
if opt.warpType=="homography":
A = tf.concat([tf.stack([X,Y,I,O,O,O,-X*(X+dX),-Y*(X+dX)],axis=-1),
tf.stack([O,O,O,X,Y,I,-X*(Y+dY),-Y*(Y+dY)],axis=-1)],1)
b = tf.expand_dims(tf.concat([X+dX,Y+dY],1),-1)
pPert = tf.matrix_solve(A,b)[:,:,0]
pPert -= tf.to_float([[1,0,0,0,1,0,0,0]])
else:
if opt.warpType=="translation":
J = np.concatenate([np.stack([I,O],axis=-1),
np.stack([O,I],axis=-1)],axis=1)
if opt.warpType=="similarity":
J = np.concatenate([np.stack([X,Y,I,O],axis=-1),
np.stack([-Y,X,O,I],axis=-1)],axis=1)
if opt.warpType=="affine":
J = np.concatenate([np.stack([X,Y,I,O,O,O],axis=-1),
np.stack([O,O,O,X,Y,I],axis=-1)],axis=1)
dXY = tf.expand_dims(tf.concat([dX,dY],1),-1)
pPert = tf.matrix_solve_ls(J,dXY)[:,:,0]
return pPert
# make training batch
def genPerturbations(opt):
with tf.name_scope("genPerturbations"):
X = np.tile(opt.canon4pts[:,0],[opt.batchSize,1])
Y = np.tile(opt.canon4pts[:,1],[opt.batchSize,1])
dX = tf.random_normal([opt.batchSize,4])*opt.pertScale \
+tf.random_normal([opt.batchSize,1])*opt.transScale
dY = tf.random_normal([opt.batchSize,4])*opt.pertScale \
+tf.random_normal([opt.batchSize,1])*opt.transScale
O = np.zeros([opt.batchSize,4],dtype=np.float32)
I = np.ones([opt.batchSize,4],dtype=np.float32)
# fit warp parameters to generated displacements
if opt.warpType=="homography":
A = tf.concat([tf.stack([X,Y,I,O,O,O,-X*(X+dX),-Y*(X+dX)],axis=-1),
tf.stack([O,O,O,X,Y,I,-X*(Y+dY),-Y*(Y+dY)],axis=-1)],1)
b = tf.expand_dims(tf.concat([X+dX,Y+dY],1),-1)
pPert = tf.matrix_solve(A,b)[:,:,0]
pPert -= tf.to_float([[1,0,0,0,1,0,0,0]])
else:
if opt.warpType=="translation":
J = np.concatenate([np.stack([I,O],axis=-1),
np.stack([O,I],axis=-1)],axis=1)
if opt.warpType=="similarity":
J = np.concatenate([np.stack([X,Y,I,O],axis=-1),
np.stack([-Y,X,O,I],axis=-1)],axis=1)
if opt.warpType=="affine":
J = np.concatenate([np.stack([X,Y,I,O,O,O],axis=-1),
np.stack([O,O,O,X,Y,I],axis=-1)],axis=1)
dXY = tf.expand_dims(tf.concat([dX,dY],1),-1)
pPert = tf.matrix_solve_ls(J,dXY)[:,:,0]
return pPert
# make training batch
def solve_linear_tf(model):
"""
Linearly solve the structure.
"""
Logger.info('Solving linear model with %d DOFs using TensorFlow...'%model.DOF)
K_bar,F_bar,index=model.K_,model.F_,model.index
K_bar=K_bar.astype(np.float32)
F_bar=F_bar.astype(np.float32)
Dvec=model.D
#Begin a new graph
if 'sess' in locals() and sess is not None:
print('Close interactive session')
sess.close()
# with tf.device('/cpu:0'):
K_init = tf.placeholder(tf.float32, shape=(model.DOF, model.DOF))
K_ = tf.Variable(K_init)
# K_ = tf.constant(K_bar,name="stiffness")
F_ = tf.Variable(np.array([F_bar.astype(np.float32)]),name="force")
# with tf.device('/cpu:0'):
D_=tf.matrix_solve(K_,tf.transpose(F_),name='displacement')
sess = tf.Session(config=tf.ConfigProto(log_device_placement=True))
print(type(K_bar),K_bar.dtype,K_bar.shape)
sess.run(tf.global_variables_initializer(),feed_dict={K_init:K_bar})
# run the op.
delta=sess.run(D_)
n_nodes=model.node_count
#fill original displacement vector
prev = 0
for idx in index:
gap=idx-prev
if gap>0:
delta=np.insert(delta,prev,[0]*gap)
prev = idx + 1
if idx==index[-1] and idx!=n_nodes-1:
delta = np.insert(delta,prev, [0]*(n_nodes*6-prev))
delta += Dvec
model.is_solved=True
return delta
