从上三角初始化对称的Theano dmatrix

Question

我正在尝试使用Theano构建一个模型，这个模型的一部分是由一个对称矩阵A来参数化的。为了确保A是对称的，我想只通过传入上三角部分的值来构建A。

相应的numpy代码可能看起来像这样：

import numpy as np

def make_symmetric(p, n):
    A = np.empty((n, n), P.dtype)
    A[np.triu_indices(n)] = p
    A.T[np.triu_indices(n)] = p

# output matrix will be (n, n)
n = 4

# parameter vector
P = np.arange(n * (n + 1) / 2)

print make_symmetric(P, n)
# [[ 0.  1.  2.  3.]
#  [ 1.  4.  5.  6.]
#  [ 2.  5.  7.  8.]
#  [ 3.  6.  8.  9.]]    

但是，由于符号张量变量不支持直接赋值，我在Theano中找不到实现这个功能的方法。

我找到的最接近的东西是theano.tensor.diag，它可以让我从对角线构建一个符号矩阵：

import theano
from theano import tensor as te

P = te.dvector('P')
D = te.diag(P)
get_D = theano.function([P], D)

print get_D(np.arange(1, 5))
# [[ 1.  0.  0.  0.]
#  [ 0.  2.  0.  0.]
#  [ 0.  0.  3.  0.]
#  [ 0.  0.  0.  4.]]

虽然还有一个theano.tensor.triu函数，但这个函数不能用来从上三角部分构建矩阵，而是返回一个数组的副本，并将下三角的元素置为零。

有没有办法从上三角部分构建一个Theano符号矩阵呢？

Answer 1

你可以使用 theano.tensor.triu 这个函数，然后把结果和它的转置相加，再减去对角线上的值。

可以直接复制粘贴的代码：

import numpy as np
import theano
import theano.tensor as T
theano.config.floatX = 'float32'

mat = T.fmatrix()
sym1 = T.triu(mat) + T.triu(mat).T
diag = T.diag(T.diagonal(mat))
sym2 = sym1 - diag

f_sym1 = theano.function([mat], sym1)
f_sym2 = theano.function([mat], sym2)

m = np.arange(9).reshape(3, 3).astype(np.float32)

print m
# [[ 0.  1.  2.]
#  [ 3.  4.  5.]
#  [ 6.  7.  8.]]
print f_sym1(m)
# [[  0.   1.   2.]
#  [  1.   8.   5.]
#  [  2.   5.  16.]]
print f_sym2(m)
# [[ 0.  1.  2.]
#  [ 1.  4.  5.]
#  [ 2.  5.  8.]]

这样做有帮助吗？这种方法需要传入一个完整的矩阵，但会忽略对角线以下的部分，并利用上三角部分来对称化。

我们还可以看看这个函数的导数。为了避免处理多维的输出，我们可以例如关注矩阵元素总和的梯度。

sum_grad = T.grad(cost=sym2.sum(), wrt=mat)
f_sum_grad = theano.function([mat], sum_grad)

print f_sum_grad(m)
# [[ 1.  2.  2.]
#  [ 0.  1.  2.]
#  [ 0.  0.  1.]]

这反映了上三角部分的元素在总和中出现了两次。

---

更新：你可以进行正常的索引操作：

n = 4
num_triu_entries = n * (n + 1) / 2

triu_index_matrix = np.zeros([n, n], dtype=int)
triu_index_matrix[np.triu_indices(n)] = np.arange(num_triu_entries)
triu_index_matrix[np.triu_indices(n)[::-1]] = np.arange(num_triu_entries)

triu_vec = T.fvector()
triu_mat = triu_vec[triu_index_matrix]

f_triu_mat = theano.function([triu_vec], triu_mat)

print f_triu_mat(np.arange(1, num_triu_entries + 1).astype(np.float32))

# [[  1.   2.   3.   4.]
#  [  2.   5.   6.   7.]
#  [  3.   6.   8.   9.]
#  [  4.   7.   9.  10.]]

---

更新：为了动态地实现这一切，一种方法是写一个符号版本的 triu_index_matrix。这可以通过一些 arange 的排列组合来完成。不过我可能在把事情搞得太复杂了。

n = T.iscalar()
n_triu_entries = (n * (n + 1)) / 2
r = T.arange(n)

tmp_mat = r[np.newaxis, :] + (n_triu_entries - n - (r * (r + 1)) / 2)[::-1, np.newaxis]
triu_index_matrix = T.triu(tmp_mat) + T.triu(tmp_mat).T - T.diag(T.diagonal(tmp_mat))

triu_vec = T.fvector()
sym_matrix = triu_vec[triu_index_matrix]

f_triu_index_matrix = theano.function([n], triu_index_matrix)
f_dynamic_sym_matrix = theano.function([triu_vec, n], sym_matrix)

print f_triu_index_matrix(5)
# [[ 0  1  2  3  4]
#  [ 1  5  6  7  8]
#  [ 2  6  9 10 11]
#  [ 3  7 10 12 13]
# [ 4  8 11 13 14]]
print f_dynamic_sym_matrix(np.arange(1., 16.).astype(np.float32), 5)
# [[  1.   2.   3.   4.   5.]
#  [  2.   6.   7.   8.   9.]
#  [  3.   7.  10.  11.  12.]
#  [  4.   8.  11.  13.  14.]
#  [  5.   9.  12.  14.  15.]]