【问题标题】:pdist for theano tensor用于 theano 张量的 pdist
【发布时间】:2014-11-11 05:48:35
【问题描述】:

我有一个theano符号矩阵

x = T.fmatrix('input')

x 稍后将由 n 暗淡 d 的向量填充(在火车时间)。

我想要pdistscipy.spatial.distance.pdist of pdist)的 theano 等价物,类似于

D = theano.pdist( x )

我怎样才能做到这一点?

x 上直接调用scipy.spatial.distance.pdist 不起作用,因为现阶段x 只是象征性的......

更新:我非常希望能够模仿 pdist“紧凑”行为:也就是说,只计算 ~1/2 nxn 条目距离矩阵。

【问题讨论】:

  • 您是在寻找 pdist 的全部通用性,还是对特定实例感兴趣,例如欧几里得距离?
  • @eickenberg 欧几里得距离将是一个不错的开始。我相信一旦我实现了这个,我就能将它推广到其他指标。

标签: python matlab scipy theano


【解决方案1】:

来自 scipy 的pdist 是不同函数的集合 - 不存在所有这些函数的 Theano 等价物。然而,每个特定的距离,作为一个封闭形式的数学表达式,都可以在 Theano 中写下来然后编译。

以minkowskip范数距离为例(复制+粘贴):

import theano
import theano.tensor as T
X = T.fmatrix('X')
Y = T.fmatrix('Y')
P = T.scalar('P')
translation_vectors = X.reshape((X.shape[0], 1, -1)) - Y.reshape((1, Y.shape[0], -1))
minkowski_distances = (abs(translation_vectors) ** P).sum(2) ** (1. / P)
f_minkowski = theano.function([X, Y, P], minkowski_distances)

注意abs调用了内置的__abs__,所以abs也是一个theano函数。我们现在可以将其与pdist 进行比较:

import numpy as np
from scipy.spatial.distance import pdist

rng = np.random.RandomState(42)
d = 20 # dimension
nX = 10
nY = 30
x = rng.randn(nX, d).astype(np.float32)
y = rng.randn(nY, d).astype(np.float32)

ps = [1., 3., 2.]

for p in ps:
    d_theano = f_minkowski(x, x, p)[np.triu_indices(nX, 1)]
    d_scipy = pdist(x, p=p, metric='minkowski')
    print "Testing p=%1.2f, discrepancy %1.3e" % (p, np.sqrt(((d_theano - d_scipy) ** 2).sum()))

这会产生

Testing p=1.00, discrepancy 1.322e-06
Testing p=3.00, discrepancy 4.277e-07
Testing p=2.00, discrepancy 4.789e-07

如您所见,对应关系是存在的,但函数f_minkowski 稍微更通用一些,因为它比较了两个可能不同的数组的行。如果将两次相同的数组作为输入传递,f_minkowski 返回一个矩阵,而pdist 返回一个没有冗余的列表。如果需要这种行为,也可以完全动态地实现,但我会在这里坚持一般情况。

但应注意一种特殊化的可能性:在 p=2 的情况下,通过二项式公式计算变得更简单,这可用于节省内存中的宝贵空间:而一般的 Minkowski 距离,如上面实现的,创建一个 3D 数组(由于避免了 for 循环和累积求和),这是禁止的,具体取决于维度 d(和 nX, nY),对于 p=2,我们可以写

squared_euclidean_distances = (X ** 2).sum(1).reshape((X.shape[0], 1)) + (Y ** 2).sum(1).reshape((1, Y.shape[0])) - 2 * X.dot(Y.T)
f_euclidean = theano.function([X, Y], T.sqrt(squared_euclidean_distances))

只使用O(nX * nY)空格而不是O(nX * nY * d)我们检查对应关系,这次是关于一般问题:

d_eucl = f_euclidean(x, y)
d_minkowski2 = f_minkowski(x, y, 2.)
print "Comparing f_minkowski, p=2 and f_euclidean: l2-discrepancy %1.3e" % ((d_eucl - d_minkowski2) ** 2).sum()

屈服

Comparing f_minkowski, p=2 and f_euclidean: l2-discrepancy 1.464e-11

【讨论】:

  • (如果这一切都在你的用例之外,请告诉我:))
  • 感谢您的详细解答。我对只获取pdist 列表的情况更感兴趣,即更紧凑的模式,其中xy 输入实际上是相同的,大约1/2 的计算没有进行。跨度>
  • 好的,应该可以提取上三角了,可能theano编译器已经理解并优化了,但我不确定。我写了一些东西来从一个上三角列表到一个完整的矩阵here,并且应该能够做相反的事情。我会检查的。
  • @eickenberg 我正在使用你的squared_euclidean_distances 实现,但它有一些负值,sqrt 包含一些NaNs。
  • 另外,请注意平方根的梯度在接近零时变得不确定。因此,我建议不要在零处剪裁,而是在稍高的值处剪裁,或者使用T.zero_grad 停止渐变。
【解决方案2】:

我以前没有使用过 Theano,但是这里有一个基于纯 Numpy 函数的解决方案(也许你将它转换为等效的 theano 函数。请注意,我在下面的表达式,所以如果 Theano 不支持它,你可能必须明确地重写它):

# X is an m-by-n matrix (rows are examples, columns are dimensions)
# D is an m-by-m symmetric matrix of pairwise Euclidean distances
a = np.sum(X**2, axis=1)
D = np.sqrt((a + a[np.newaxis].T) - 2*np.dot(X, X.T))

它基于以下事实:||u-v||^2 = ||u||^2 + ||v||^2 - 2*u.v。 (我使用 MATLAB 在我的 previous answers 中展示了这个)

这是与 Scipy 现有函数的比较:

import numpy as np
from scipy.spatial.distance import pdist, squareform

def my_pdist(X):
    a = np.sum(X**2, axis=1)
    D = np.sqrt((a + a[np.newaxis].T) - 2*np.dot(X, X.T))
    return D

def scipy_pdist(X):
    D = squareform(pdist(X, metric='euclidean'))
    return D    

X = np.random.rand(5, 3)
D1 = my_pdist(X)
D2 = scipy_pdist(X)

差异应该可以忽略不计,接近机器 epsilon (np.spacing(1)):

>>> np.linalg.norm(D1-D2)
8.5368137554718277e-16

HTH


编辑:

这是另一个使用单个循环的实现:

def my_pdist_compact(X):
    D = np.empty(shape=[0,0], dtype=X.dtype)
    for i in range(X.shape[0]-1):
        D = np.append(D, np.sqrt(np.sum((X[i,] - X[i+1:,])**2, axis=1)))
    return D

有点等价的 MATLAB 代码:

function D = my_pdist_compact(X)
    n = size(X,1);
    D = cell(n-1,1);
    for i=1:n-1
        D{i} = sqrt(sum(bsxfun(@minus, X(i,:), X(i+1:end,:)).^2, 2));
    end
    D = vertcat(D{:});
end

这将返回紧凑形式的成对距离(对称矩阵的上三角部分)。这与pdist 的输出相同。使用squareform 将其转换为完整矩阵。

>>> d1 = my_pdist_compact(X)
>>> d2 = pdist(X)    # from scipy.spatial.distance
>>> (d1 == d2).all()
True

我会留给你看看是否可以使用 Theano 编写等效的loop(请参阅theano.scan)!

【讨论】:

  • 谢谢!但是,我正在寻找一种方法来保存 squareform 结果并直接获得紧凑的结果。
  • 此外,我正在寻找一种更“theano”的方式来执行此操作,因为编译器可能更容易优化...
  • 如果这个错误忽略我(就像我说我从未使用过 Theano),但你不能只使用常规循环实现 pdist (你可以迭代行甚至列)。据我了解,代码被编译成 C 代码,所以解释语言的循环开销就没有了。如果它足够聪明,它应该能够从 BLAS 调用向量化函数。
  • 恐怕我的theano技能对于构建符号嵌套循环来说太基础了......
  • Theano 在这方面更像 SymPy:它为您提供了一个它知道如何解释的函数目录,并且您可以象征性地使用通过该目录的组合可以访问的各种表达式。但它不擅长吸收目录之外的东西。使用 Theano 时的一个好的经验法则是:问问自己,“我希望 SymPy 为我处理这种计算吗?”如果答案是否定的,那么您可能(并非总是)不恰当地使用 Theano。这并不是说你不能让它工作。只是有更合适的方法。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2016-11-06
  • 1970-01-01
  • 2015-01-19
  • 2015-11-14
  • 1970-01-01
相关资源
最近更新 更多