【问题标题】:Linear projection matrix by solving the procrustes problem通过解决 procrustes 问题的线性投影矩阵
【发布时间】:2019-09-04 13:15:35
【问题描述】:

我有两个矩阵:

target = np.array([[1, 1, 1, 1, 1],
               [2, 2, 2, 2, 2],
               [3, 3, 3, 3, 3]])

source = np.array([[11, 11, 11, 11, 11],
               [22, 22, 22, 22, 22],
               [33, 33, 33, 33, 33]])

我想创建一个转换矩阵,将矩阵投影到目标矩阵。

我发现 Scipy 库提供了一个功能来做到这一点:

from scipy.spatial import procrustes
mtx1, mtx2, disparity = procrustes(target, source)

基于documentation,它说:

因此,mtx2 是投影矩阵。

如果我有另一个数据,我想使用 Scipy 用于将 矩阵投影到目标一个?

我如何使用 Scipy 做到这一点?

【问题讨论】:

    标签: python matrix scipy procrustes


    【解决方案1】:

    您需要修改函数才能返回转换矩阵(R)。

    移除 cmets 后的源代码如下所示:

    def procrustes(data1, data2):
        mtx1 = np.array(data1, dtype=np.double, copy=True)
        mtx2 = np.array(data2, dtype=np.double, copy=True)
    
        if mtx1.ndim != 2 or mtx2.ndim != 2:
            raise ValueError("Input matrices must be two-dimensional")
        if mtx1.shape != mtx2.shape:
            raise ValueError("Input matrices must be of same shape")
        if mtx1.size == 0:
            raise ValueError("Input matrices must be >0 rows and >0 cols")
    
        # translate all the data to the origin
        mtx1 -= np.mean(mtx1, 0)
        mtx2 -= np.mean(mtx2, 0)
    
        norm1 = np.linalg.norm(mtx1)
        norm2 = np.linalg.norm(mtx2)
    
        if norm1 == 0 or norm2 == 0:
            raise ValueError("Input matrices must contain >1 unique points")
    
        # change scaling of data (in rows) such that trace(mtx*mtx') = 1
        mtx1 /= norm1
        mtx2 /= norm2
    
        # transform mtx2 to minimize disparity
        R, s = orthogonal_procrustes(mtx1, mtx2)
        mtx2 = np.dot(mtx2, R.T) * s    # HERE, the projected mtx2 is estimated.
    
        # measure the dissimilarity between the two datasets
        disparity = np.sum(np.square(mtx1 - mtx2))
    
        return mtx1, mtx2, disparity, R
    

    来源:https://github.com/scipy/scipy/blob/v1.3.0/scipy/spatial/_procrustes.py#L17-L132

    【讨论】:

    • Thanx @serafeim ,这就是我要找的。最后一个问题,如何将 mtx2 矩阵恢复为源形状?
    • mtx2_proj = np.dot(mtx2, R.T) 所以mtx2 = np.dot(mtx2_proj, R)
    • 不加 /s?我问你是因为我注意到当我像你一样这样做时,我得到了与原始“相似”但在零附近标准化(因为代码中的第一步,标准化)。
    • 对不起,是的。您需要在代码中重命名mtx2,因为在应用规范化时它会被覆盖。
    猜你喜欢
    • 1970-01-01
    • 2011-08-31
    • 1970-01-01
    • 2016-02-15
    • 2020-10-22
    • 2012-03-27
    • 1970-01-01
    • 2015-11-01
    • 2013-08-03
    相关资源
    最近更新 更多