如何计算成对距离函数的雅可比行列式（`scipy.spatial.pdist`）

Question

语境

我是netgraph的作者和维护者，netgraph 是一个用于创建网络可视化的 Python 库。N我目前正在尝试优化一个例程，该例程计算每条边都有定义长度的网络的一组节点位置。可以在此处找到示例。

问题

该例程的核心是scipy.optimize.minimize计算使节点之间的总距离最大化的位置：

def cost_function(positions):
    return 1. / np.sum((pdist(positions.reshape((-1, 2))))**power)

result = minimize(cost_function, initial_positions.flatten(), method='SLSQP',
                  jac="2-point", constraints=[nonlinear_constraint])

• positions是 (x, y) 元组的（已解开的）numpy 数组。
• power是一个较小的数字，限制了大距离的影响（以鼓励紧凑的节点布局），但出于本问题的目的，可以假设为 1。
• pdist是 中的成对距离函数scipy.spatial。

最小化（/最大化）使用以下非线性约束进行约束：

lower_bounds = ... # (squareform of an) (N, N) distance matrix of the sum of node sizes (i.e. nodes should not overlap)
upper_bounds = ... # (squareform of an) (N, N) distance matrix constructed from the given edge lengths

def constraint_function(positions):
    positions = np.reshape(positions, (-1, 2))
    return pdist(positions)

nonlinear_constraint = NonlinearConstraint(constraint_function, lb=lower_bounds, ub=upper_bounds, jac='2-point')

对于玩具示例，优化可以正确且快速地完成。然而，即使对于小型网络，运行时间也相当糟糕。我当前的实现使用有限差分来近似梯度 ( jac='2-point')。为了加快计算速度，我想显式计算雅可比行列式。

在几篇 Math Stackexchange 帖子（1、2 ）之后，我计算了成对距离函数的雅可比行列式，如下所示：

    def delta_constraint(positions):
        positions = np.reshape(positions, (-1, 2))
        total_positions = positions.shape[0]
        delta = positions[np.newaxis, :, :] - positions[:, np.newaxis, :]
        distance = np.sqrt(np.sum(delta ** 2, axis=-1))
        jac = delta / distance[:, :, np.newaxis]
        squareform_indices = np.triu_indices(total_positions, 1)
        return jac[squareform_indices]

nonlinear_constraint = NonlinearConstraint(constraint_function, lb=lower_bounds, ub=upper_bounds, jac=delta_constraint)

但是，这会导致ValueError，因为输出的形状不正确。对于三角形示例，预期输出形状为 (3, 6)，而上面的函数返回 (3, 2) 数组（即 3 个成对距离乘以 2 维）。对于正方形，预期输出为 (6, 8)，而实际输出为 (6, 2)。任何有助于为和 的jac参数导出实现正确的可调用对象的帮助将不胜感激。NonlinearConstraintminimize

笔记

我想避免使用 autograd/jax/numdifftools （如在这个问题中），因为我想保持我的库的依赖项数量较小。

最小工作示例

#!/usr/bin/env python
"""
Create a node layout with fixed edge lengths but unknown node positions.
"""

import numpy as np

from scipy.optimize import minimize, NonlinearConstraint
from scipy.spatial.distance import pdist, squareform


def get_geometric_node_layout(edges, edge_length, node_size=0., power=0.2, maximum_iterations=200, origin=(0, 0), scale=(1, 1)):
    """Node layout for defined edge lengths but unknown node positions.

    Node positions are determined through non-linear optimisation: the
    total distance between nodes is maximised subject to the constraint
    imposed by the edge lengths, which are used as upper bounds.
    If provided, node sizes are used to set lower bounds.

    Parameters
    ----------
    edges : list
        The edges of the graph, with each edge being represented by a (source node ID, target node ID) tuple.
    edge_lengths : dict
        Mapping of edges to their lengths.
    node_size : scalar or dict, default 0.
        Size (radius) of nodes.
        Providing the correct node size minimises the overlap of nodes in the graph,
        which can otherwise occur if there are many nodes, or if the nodes differ considerably in size.
    power : float, default 0.2.
        The cost being minimised is the inverse of the sum of distances.
        The power parameter is the exponent applied to each distance before summation.
        Large values result in positions that are stretched along one axis.
        Small values decrease the influence of long distances on the cost
        and promote a more compact layout.
    maximum_iterations : int
        Maximum number of iterations of the minimisation.
    origin : tuple, default (0, 0)
        The (float x, float y) coordinates corresponding to the lower left hand corner of the bounding box specifying the extent of the canvas.
    scale : tuple, default (1, 1)
        The (float x, float y) dimensions representing the width and height of the bounding box specifying the extent of the canvas.

    Returns
    -------
    node_positions : dict
        Dictionary mapping each node ID to (float x, float y) tuple, the node position.

    """
    # TODO: assert triangle inequality

    # TODO: assert that the edges fit within the canvas dimensions

    # ensure that graph is bi-directional
    edges = edges + [(target, source) for (source, target) in edges] # forces copy
    edges = list(set(edges))

    # upper bound: pairwise distance matrix with unknown distances set to the maximum possible distance given the canvas dimensions
    lengths = []
    for (source, target) in edges:
        if (source, target) in edge_length:
            lengths.append(edge_length[(source, target)])
        else:
            lengths.append(edge_length[(target, source)])

    sources, targets = zip(*edges)
    nodes = sources + targets
    unique_nodes = set(nodes)
    indices = range(len(unique_nodes))
    node_to_idx = dict(zip(unique_nodes, indices))
    source_indices = [node_to_idx[source] for source in sources]
    target_indices = [node_to_idx[target] for target in targets]

    total_nodes = len(unique_nodes)
    max_distance = np.sqrt(scale[0]**2 + scale[1]**2)
    distance_matrix = np.full((total_nodes, total_nodes), max_distance)
    distance_matrix[source_indices, target_indices] = lengths
    distance_matrix[np.diag_indices(total_nodes)] = 0
    upper_bounds = squareform(distance_matrix)

    # lower bound: sum of node sizes
    if isinstance(node_size, (int, float)):
        sizes = node_size * np.ones((total_nodes))
    elif isinstance(node_size, dict):
        sizes = np.array([node_size[node] if node in node_size else 0. for node in unique_nodes])

    sum_of_node_sizes = sizes[np.newaxis, :] + sizes[:, np.newaxis]
    sum_of_node_sizes -= np.diag(np.diag(sum_of_node_sizes)) # squareform requires zeros on diagonal
    lower_bounds = squareform(sum_of_node_sizes)

    def cost_function(positions):
        return 1. / np.sum((pdist(positions.reshape((-1, 2))))**power)

    def constraint_function(positions):
        positions = np.reshape(positions, (-1, 2))
        return pdist(positions)

    initial_positions = _initialise_geometric_node_layout(edges)
    nonlinear_constraint = NonlinearConstraint(constraint_function, lb=lower_bounds, ub=upper_bounds, jac='2-point')
    result = minimize(cost_function, initial_positions.flatten(), method='SLSQP',
                      jac="2-point", constraints=[nonlinear_constraint], options=dict(maxiter=maximum_iterations))

    if not result.success:
        print("Warning: could not compute valid node positions for the given edge lengths.")
        print(f"scipy.optimize.minimize: {result.message}.")

    node_positions_as_array = result.x.reshape((-1, 2))
    node_positions = dict(zip(unique_nodes, node_positions_as_array))
    return node_positions


def _initialise_geometric_node_layout(edges):
    sources, targets = zip(*edges)
    total_nodes = len(set(sources + targets))
    return np.random.rand(total_nodes, 2)


if __name__ == '__main__':

    import matplotlib.pyplot as plt

    def plot_graph(edges, node_layout):
        # poor man's graph plotting
        fig, ax = plt.subplots()
        for source, target in edges:
            x1, y1 = node_layout[source]
            x2, y2 = node_layout[target]
            ax.plot([x1, x2], [y1, y2], color='darkgray')
        ax.set_aspect('equal')

    ################################################################################
    # triangle with right angle

    edges = [
        (0, 1),
        (1, 2),
        (2, 0)
    ]

    lengths = {
        (0, 1) : 3,
        (1, 2) : 4,
        (2, 0) : 5,
    }

    pos = get_geometric_node_layout(edges, lengths, node_size=0)

    plot_graph(edges, node_layout=pos)

    plt.show()

    ################################################################################
    # square

    edges = [
        (0, 1),
        (1, 2),
        (2, 3),
        (3, 0),
    ]

    lengths = {
        (0, 1) : 0.5,
        (1, 2) : 0.5,
        (2, 3) : 0.5,
        (3, 0) : 0.5,
    }

    pos = get_geometric_node_layout(edges, lengths, node_size=0)

    plot_graph(edges, node_layout=pos)

    plt.show()

编辑：计时的现实用例

下面是一个更实际的用例，我用它来计时我的代码。我已经将 @adrianop01 的雅可比行列式计算纳入约束中。它还包括高级初始化。它需要额外的依赖项networkx和netgraph，这两个依赖项都可以通过 pip 安装。

#!/usr/bin/env python
"""
Create a node layout with fixed edge lengths but unknown node positions.
"""

import numpy as np

from itertools import combinations
from scipy.optimize import minimize, NonlinearConstraint
from scipy.spatial.distance import pdist, squareform

from netgraph._node_layout import _rescale_to_frame


def get_geometric_node_layout(edges, edge_length, node_size=0., power=0.2, maximum_iterations=200, origin=(0, 0), scale=(1, 1)):
    """Node layout for defined edge lengths but unknown node positions.

    Node positions are determined through non-linear optimisation: the
    total distance between nodes is maximised subject to the constraint
    imposed by the edge lengths, which are used as upper bounds.
    If provided, node sizes are used to set lower bounds.

    Parameters
    ----------
    edges : list
        The edges of the graph, with each edge being represented by a (source node ID, target node ID) tuple.
    edge_lengths : dict
        Mapping of edges to their lengths.
    node_size : scalar or dict, default 0.
        Size (radius) of nodes.
        Providing the correct node size minimises the overlap of nodes in the graph,
        which can otherwise occur if there are many nodes, or if the nodes differ considerably in size.
    power : float, default 0.2.
        The cost being minimised is the inverse of the sum of distances.
        The power parameter is the exponent applied to each distance before summation.
        Large values result in positions that are stretched along one axis.
        Small values decrease the influence of long distances on the cost
        and promote a more compact layout.
    maximum_iterations : int
        Maximum number of iterations of the minimisation.
    origin : tuple, default (0, 0)
        The (float x, float y) coordinates corresponding to the lower left hand corner of the bounding box specifying the extent of the canvas.
    scale : tuple, default (1, 1)
        The (float x, float y) dimensions representing the width and height of the bounding box specifying the extent of the canvas.

    Returns
    -------
    node_positions : dict
        Dictionary mapping each node ID to (float x, float y) tuple, the node position.

    """
    # TODO: assert triangle inequality

    # TODO: assert that the edges fit within the canvas dimensions

    # ensure that graph is bi-directional
    edges = edges + [(target, source) for (source, target) in edges] # forces copy
    edges = list(set(edges))

    # upper bound: pairwise distance matrix with unknown distances set to the maximum possible distance given the canvas dimensions
    lengths = []
    for (source, target) in edges:
        if (source, target) in edge_length:
            lengths.append(edge_length[(source, target)])
        else:
            lengths.append(edge_length[(target, source)])

    sources, targets = zip(*edges)
    nodes = sources + targets
    unique_nodes = set(nodes)
    indices = range(len(unique_nodes))
    node_to_idx = dict(zip(unique_nodes, indices))
    source_indices = [node_to_idx[source] for source in sources]
    target_indices = [node_to_idx[target] for target in targets]

    total_nodes = len(unique_nodes)
    max_distance = np.sqrt(scale[0]**2 + scale[1]**2)
    distance_matrix = np.full((total_nodes, total_nodes), max_distance)
    distance_matrix[source_indices, target_indices] = lengths
    distance_matrix[np.diag_indices(total_nodes)] = 0
    upper_bounds = squareform(distance_matrix)

    # lower bound: sum of node sizes
    if isinstance(node_size, (int, float)):
        sizes = node_size * np.ones((total_nodes))
    elif isinstance(node_size, dict):
        sizes = np.array([node_size[node] if node in node_size else 0. for node in unique_nodes])

    sum_of_node_sizes = sizes[np.newaxis, :] + sizes[:, np.newaxis]
    sum_of_node_sizes -= np.diag(np.diag(sum_of_node_sizes)) # squareform requires zeros on diagonal
    lower_bounds = squareform(sum_of_node_sizes)
    invalid = lower_bounds > upper_bounds
    lower_bounds[invalid] = upper_bounds[invalid] - 1e-8

    def cost_function(positions):
        # return -np.sum((pdist(positions.reshape((-1, 2))))**power)
        return 1. / np.sum((pdist(positions.reshape((-1, 2))))**power)

    def cost_jacobian(positions):
        # TODO
        pass

    def constraint_function(positions):
        positions = np.reshape(positions, (-1, 2))
        return pdist(positions)

    # adapted from /sf/answers/5260807681/
    total_pairs = int((total_nodes - 1) * total_nodes / 2)
    source_indices, target_indices = np.array(list(combinations(range(total_nodes), 2))).T # node order thus (0,1) ... (0,N-1), (1,2),...(1,N-1),...,(N-2,N-1)
    rows = np.repeat(np.arange(total_pairs).reshape(-1, 1), 2, axis=1)
    source_columns = np.vstack((source_indices*2, source_indices*2+1)).T
    target_columns = np.vstack((target_indices*2, target_indices*2+1)).T

    def constraint_jacobian(positions):
        positions = np.reshape(positions, (-1, 2))
        pairwise_distances = constraint_function(positions)
        jac = np.zeros((total_pairs, 2 * total_nodes))
        jac[rows, source_columns] = (positions[source_indices] - positions[target_indices]) / pairwise_distances.reshape((-1, 1))
        jac[rows, target_columns] = -jac[rows, source_columns]
        return jac

    initial_positions = _initialise_geometric_node_layout(edges, edge_length)
    nonlinear_constraint = NonlinearConstraint(constraint_function, lb=lower_bounds, ub=upper_bounds, jac=constraint_jacobian)
    result = minimize(cost_function, initial_positions.flatten(), method='SLSQP',
                      jac='2-point', constraints=[nonlinear_constraint], options=dict(maxiter=maximum_iterations))
    # result = minimize(cost_function, initial_positions.flatten(), method='trust-constr',
    #                   jac=cost_jacobian, constraints=[nonlinear_constraint])

    if not result.success:
        print("Warning: could not compute valid node positions for the given edge lengths.")
        print(f"scipy.optimize.minimize: {result.message}.")

    node_positions_as_array = result.x.reshape((-1, 2))
    node_positions_as_array = _rescale_to_frame(node_positions_as_array, np.array(origin), np.array(scale))
    node_positions = dict(zip(unique_nodes, node_positions_as_array))
    return node_positions


# # slow
# def _initialise_geometric_node_layout(edges, edge_length=None):
#     sources, targets = zip(*edges)
#     total_nodes = len(set(sources + targets))
#     return np.random.rand(total_nodes, 2)

# much faster
def _initialise_geometric_node_layout(edges, edge_length=None):
    """Initialises the node positions using the FR algorithm with weights.
    Shorter edges are given a larger weight such that the nodes experience a strong attractive force."""

    from netgraph import get_fruchterman_reingold_layout
    if edge_length:
        edge_weight = dict()
        for edge, length in edge_length.items():
            edge_weight[edge] = 1 / length
    else:
        edge_weight = None
    node_positions = get_fruchterman_reingold_layout(edges)
    return np.array(list(node_positions.values()))


if __name__ == '__main__':

    from time import time
    import matplotlib.pyplot as plt
    import networkx as nx # pip install networkx

    from netgraph import Graph # pip install netgraph


    fig, (ax1, ax2) = plt.subplots(1, 2)

    g = nx.random_geometric_graph(50, 0.3, seed=2)
    node_positions = nx.get_node_attributes(g, 'pos')
    plot_instance = Graph(g,
                          node_layout=node_positions,
                          node_size=1, # netgraph rescales node sizes by 0.01
                          node_edge_width=0.1,
                          edge_width=0.1,
                          ax=ax1,
    )
    ax1.axis([0, 1, 0, 1])
    ax1.set_title('Original node positions')

    def get_euclidean_distance(p1, p2):
        return np.sqrt(np.sum((np.array(p1)-np.array(p2))**2))

    edge_length = dict()
    for (source, target) in g.edges:
        edge_length[(source, target)] = get_euclidean_distance(node_positions[source], node_positions[target])

    tic = time()
    new_node_positions = get_geometric_node_layout(list(g.edges), edge_length, node_size=0.01)
    toc = time()

    print(f"Time elapsed : {toc-tic}")

    Graph(g,
          node_layout=new_node_positions,
          node_size=1,
          node_edge_width=0.1,
          edge_width=0.1,
          ax=ax2,
    )
    ax2.axis([0, 1, 0, 1])
    ax2.set_title('Reconstructed node positions')

    plt.show()

---

第二次编辑

以下是我在测试 @spinkus 和相关解决方案时获得的一些初步结果。我对他的代码的实现如下所示：

def cost_function(positions):
    return -np.sum((pdist(positions.reshape((-1, 2))))**2)

def cost_jacobian(positions):
    positions = positions.reshape(-1, 2)
    delta = positions[np.newaxis, :] - positions[:, np.newaxis]
    jac = -2 * np.sum(delta, axis=0)
    return jac.ravel()

不幸的是，这个成本函数需要更长的时间才能收敛：5 次最好的情况下需要 13 秒，并且时间差异很大（最多一分钟）。这与我使用显式雅可比行列式还是使用有限差分方法近似它无关。此外，最小化通常会因“scipy.optimize.minimize：不等式约束不兼容”而提前结束。和“scipy.optimize.minimize：线性搜索的正向导数”。我的赌注（尽管我没有什么证据支持）是成本的绝对值很重要。我原来的成本函数在价值和绝对值上都下降了，而最小化增加了@spinkus成本函数的绝对值（但是，请参阅下面@spinkus的优秀评论，为什么这可能有点转移注意力并导致不太准确的解决方案）。

我也理解（我认为）为什么我原来的成本函数不适合计算雅可比行列式。设power为 0.5，则成本函数和雅可比行列式采用以下形式（除非我的代数又错了）：

def cost_function(positions):
    return 1. / np.sum((pdist(positions.reshape((-1, 2))))**0.5)

def cost_jacobian(positions):
    positions = positions.reshape(-1, 2)
    delta = positions[np.newaxis, :] - positions[:, np.newaxis]
    distance = np.sqrt(np.sum(delta**2, axis=-1))
    denominator = -2 * np.sqrt(delta) * distance[:, :, np.newaxis]
    denominator[np.diag_indices_from(denominator[:, :, 0]),:] = 1
    jac = 1 / denominator
    return np.sum(jac, axis=0).ravel() - 1

有问题的术语是sqrt(delta)，其中delta是所有点之间的向量。忽略对角线，该矩阵中的一半条目必然为负，因此无法计算雅可比行列式。

然而，该电源的目的只是为了降低长距离对成本的重要性。任何具有递减导数的单调递增函数都可以。使用log(x + 1)电源代替电源会产生以下功能：

def cost_function(positions):
    return 1 / np.sum(np.log(pdist(positions.reshape((-1, 2))) + 1))

def cost_jacobian(positions):
    positions = positions.reshape(-1, 2)
    delta = positions[np.newaxis, :] - positions[:, np.newaxis]
    distance2 = np.sum(delta**2, axis=-1)
    distance2[np.diag_indices_from(distance2)] = 1
    jac = -delta / (distance2 + np.sqrt(distance2))[..., np.newaxis]
    return np.sum(jac, axis=0).ravel()

使用有限差分近似，最小化会在 0.5 秒内终止。然而，对于显式雅可比行列式，最佳运行时间为 4 秒，尽管仍然存在很大的差异，运行时间增加了一分钟甚至更长。

太尔；博士。

我仍然不明白为什么最小化不能使用显式雅可比行列式运行得更快。

Answer 1

此实现根据与 OP 的讨论，计算所有点对的约束函数的雅可比行列式。np 数组向量化代码可能并不完美，因此我欢迎进一步评论基于雅可比公式的代码细化。

雅可比矩阵的推导

2d 点可以定义为向量：

雅可比矩阵为（M 行，N 列，M 是唯一 2 点对的数量，N 是唯一点的数量）：

对于雅可比矩阵的每个单独元素，我们有以下三种情况：

1. 具体点坐标x/y是2-范数中的点1 x/y。
2. 具体点坐标x/y是2-范数中的点2x/y。
3. 具体点坐标x/y与当前的2-范数无关。

因此，我们期望雅可比行列式是一个包含大量零的稀疏矩阵，每行最多有 4 个非零项。

代码

该代码是不言自明的。我们将雅可比矩阵作为 MxNnp.zeros矩阵开始，并且仅更新与当前 2 范数函数/点对相关的那些条目（因此，每行 4 次更新）。

from itertools import combinations

n_pt = len(initial_positions) #N
n_ptpair = len(upper_bounds) #total number of pointpairs, M

idx_pts= np.array(list(combinations(range(n_pt),2))) #point id order thus in (0,1) ... (0,N-1), (1,2),...(1,N-1),...,(N-2,N-1)
idx_pt1= np.array(idx_pts[:,0]) 
idx_pt2= np.array(idx_pts[:,1])

row_idx = np.repeat(np.arange(n_ptpair).reshape(-1,1),2,axis=1)
col1_idx = np.vstack((idx_pt1*2,idx_pt1*2+1)).T
col2_idx = np.vstack((idx_pt2*2,idx_pt2*2+1)).T

def delta_constraint(positions):
  positions = np.reshape(positions, (-1, 2))
  pairdist = constraint_function(positions) #pairwise R2 distance between each point pair
  jac = np.zeros((n_ptpair,2*n_pt)) #(M,(x0,y0,x1,y1,...,xc,yc...,xN,yN))
  jac[row_idx,col1_idx] = (positions[idx_pt1]-positions[idx_pt2])/pairdist.reshape((-1,1))
  jac[row_idx,col2_idx] = -jac[row_idx,col1_idx]
  return jac

结果

1.代码中的三角形

2.代码中的Square

3.比较复杂的图

  edges = [
        (0, 1),
        (1, 2),
        (2, 3),
        (3, 0),
        (3, 1),
        (4, 1),
        (5, 1),
        (5, 2),
    ]

    lengths = {
        (0, 1) : 0.5,
        (1, 2) : 0.5,
        (2, 3) : 0.5,
        (3, 0) : 0.5,
        (3, 1) : 0.8,
        (4, 1) : 0.8,
        (5, 1) : 0.2,
        (5, 2) : 1,
    }