启发式选择使点积最大化的五列数组

Posted 2023-03-12

【中文标题】启发式选择使点积最大化的五列数组

【英文标题】：Heuristic to choose five column arrays that maximise the dot product

【发布时间】：2021-10-13 03:23:09

【问题描述】：

我有一个稀疏的 60000x10000 矩阵 M，其中每个元素都是 1 或 0。矩阵中的每一列都是不同的信号组合（即 1 和 0）。我想从 M 中选择五个列向量并取它们的 Hadamard（即元素）乘积；我将结果向量称为策略向量。在这一步之后，我计算此策略向量与目标向量（不变）的点积。目标向量用 1 和 -1 填充，这样策略向量的特定行中为 1 会得到奖励或惩罚。

我是否可以使用一些启发式或线性代数方法来帮助我从矩阵 M 中选择五个向量，从而产生高点积？我没有任何使用 Google 的经验OR工具也不是Scipy的优化方法，所以我不太确定它们是否可以应用于我的问题。对此的建议将不胜感激！ :)

注意：作为解给出的五个列向量不需要是最优的；我宁愿有一些不需要几个月/几年才能运行的东西。

【问题讨论】：

有趣的问题。五个任意列向量还是五个相邻列向量？也许您可以展示一个非常小的示例来说明您想要实现的目标？

来自 M c1、c2、c3、c4 和 c5 的五个任意列。每列用 1 和 0 填充，因此 1 是信号。然后计算这些列的元素运算，因为这类似于“和”运算符：所有列必须在一行中有一个 1，才能在该行的结果向量 s 中有一个 1。我想最大化 s 和目标向量 s* 的点积。计算可能导致策略 s 的所有可能的列组合是不可行的，所以我想要一个更快的解决方案（不需要是绝对最优的解决方案）。如果有什么不清楚的地方请告诉我:)

从 10,000 个元素中选择 5 个元素可以得到大约 10^18 种可能的组合。除非您能以某种方式将其转化为凸问题，否则这种规模的优化选项并不多。

【参考方案1】：

首先，感谢您提出一个好问题。我不会经常练习 numpy。另外，我没有太多在 SE 上发帖的经验，因此欢迎任何与答案相关的反馈、代码批评和意见。

这最初是为了寻找最佳解决方案的尝试，但我没有设法处理复杂性。但是，该算法应该为您提供一个可能被证明足够的贪心解决方案。

Colab Notebook (Python code + Octave validation)

核心理念

注意：在运行时，我已经转置了矩阵。所以，问题中的列向量对应算法中的行向量。

请注意，您可以一次将目标与一个向量相乘，从而有效地得到一个新目标，但其中包含一些 0s。这些永远不会改变，因此您可以通过在进一步的计算中完全删除那些行（算法中的列）来过滤掉一些计算 - 无论是从目标还是矩阵中。 - 然后你又得到了一个有效的目标（只有1s 和-1）。

这就是算法的基本思想。给定：

n: 你需要选择的向量数量 b: 要检查的最佳向量的数量 m: 检查一个向量的矩阵运算的复杂度

执行指数复杂的O((n*m)^b) 深度优先搜索，但通过减小目标/矩阵大小来降低更深层计算的复杂性，同时通过一些启发式方法减少一些搜索路径。

使用的启发式方法

目前为止取得的最好成绩在每个递归步骤中都是已知的。计算一个乐观向量（将-1 转为0）并检查仍然可以达到什么分数。不要在分数不能超过的关卡中搜索。

如果矩阵中的最佳向量具有均匀分布的1s 和0s，这将毫无用处。乐观的分数太高了。但是，稀疏度越高，效果越好。

忽略重复。基本上，不要检查同一层中的重复向量。因为我们减少了矩阵大小，所以在更深的递归级别中出现重复的机会会增加。

关于启发式的进一步思考 最有价值的是那些在一开始就消除了向量选择的那些。就它们对目标的影响而言，可能有一种方法可以找到比其他向量更差或相等的向量。比如说，如果v1 与v2 仅相差一个额外的1，并且目标在该行中有一个-1，那么v1 比v2 差或等于。

问题是我们需要找到超过 1 个向量，而不能轻易丢弃其余的。如果我们有 10 个向量，每个向量都比前一个差或相等，我们仍然必须在开始时保留 5 个（以防它们仍然是最佳选择），然后在下一个递归级别保留 4 个，在下面的递归级别保留 3 个等。

也许可以生成一棵树并将其传递给递归？不过，这无助于在一开始就缩小搜索空间......也许只考虑更差或相等链中的 1 或 2 个向量会有所帮助？这将探索更多样化的解决方案，但并不能保证它会更优化。

警告：请注意，示例中的 MATRIX 和 TARGET 在 int8 中。如果将这些用于点积，它将溢出。虽然我认为算法中的所有操作都在创建新变量，所以不受影响。

代码

# Given:
TARGET = np.random.choice([1, -1], size=60000).astype(np.int8)
MATRIX = np.random.randint(0, 2, size=(10000,60000), dtype=np.int8)

# Tunable - increase to search more vectors, at the cost of time.
# Performs better if the best vectors in the matrix are sparse
MAX_BRANCHES = 3   # can give more for sparser matrices

# Usage
score, picked_vectors_idx = pick_vectors(TARGET, MATRIX, 5)

# Function
def pick_vectors(init_target, init_matrix, vectors_left_to_pick: int, best_prev_result=float("-inf")):
    assert vectors_left_to_pick >= 1
    if init_target.shape == (0, ) or len(init_matrix.shape) <= 1 or init_matrix.shape[0] == 0 or init_matrix.shape[1] == 0:
        return float("inf"), None
    target = init_target.copy()
    matrix = init_matrix.copy()

    neg_matrix = np.multiply(target, matrix)
    neg_matrix_sum = neg_matrix.sum(axis=1)

    if vectors_left_to_pick == 1:
        picked_id = np.argmax(neg_matrix_sum)
        score = neg_matrix[picked_id].sum()
        return score, [picked_id]

    else:
        sort_order = np.argsort(neg_matrix_sum)[::-1]
        sorted_sums = neg_matrix_sum[sort_order]
        sorted_neg_matrix = neg_matrix[sort_order]
        sorted_matrix = matrix[sort_order]

        best_score = best_prev_result
        best_picked_vector_idx = None

        # Heuristic 1 (H1) - optimistic target.
        # Set a maximum score that can still be achieved
        optimistic_target = target.copy()
        optimistic_target[target == -1] = 0
        if optimistic_target.sum() <= best_score:
            # This check can be removed - the scores are too high at this point
            return float("-inf"), None

        # Heuristic 2 (H2) - ignore duplicates
        vecs_tried = set()

        # MAIN GOAL:   for picked_id, picked_vector in enumerate(sorted_matrix):
        for picked_id, picked_vector in enumerate(sorted_matrix[:MAX_BRANCHES]):
            # H2
            picked_tuple = tuple(picked_vector)
            if picked_tuple in vecs_tried:
                continue
            else:
                vecs_tried.add(picked_tuple)

            # Discard picked vector
            new_matrix = np.delete(sorted_matrix, picked_id, axis=0)
            
            # Discard matrix and target rows where vector is 0
            ones = np.argwhere(picked_vector == 1).squeeze()
            new_matrix = new_matrix[:, ones]
            new_target = target[ones]
            if len(new_matrix.shape) <= 1 or new_matrix.shape[0] == 0:
                return float("-inf"), None

            # H1: Do not compute if best score cannot be improved
            new_optimistic_target = optimistic_target[ones]
            optimistic_matrix = np.multiply(new_matrix, new_optimistic_target)
            optimistic_sums = optimistic_matrix.sum(axis=1)
            optimistic_viable_vector_idx = optimistic_sums > best_score
            if optimistic_sums.max() <= best_score:
                continue
            new_matrix = new_matrix[optimistic_viable_vector_idx]
         
            score, next_picked_vector_idx = pick_vectors(new_target, new_matrix, vectors_left_to_pick - 1, best_prev_result=best_score)
            
            if score <= best_score:
                continue

            # Convert idx of trimmed-down matrix into sorted matrix IDs
            for i, returned_id in enumerate(next_picked_vector_idx):
                # H1: Loop until you hit the required number of 'True'
                values_passed = 0
                j = 0
                while True:
                    value_picked: bool = optimistic_viable_vector_idx[j]
                    if value_picked:
                        values_passed += 1
                        if values_passed-1 == returned_id:
                            next_picked_vector_idx[i] = j
                            break
                    j += 1

                # picked_vector index
                if returned_id >= picked_id:
                    next_picked_vector_idx[i] += 1

            best_score = score

            # Convert from sorted matrix to input matrix IDs before returning
            matrix_id = sort_order[picked_id]
            next_picked_vector_idx = [sort_order[x] for x in next_picked_vector_idx]
            best_picked_vector_idx = [matrix_id] + next_picked_vector_idx

        return best_score, best_picked_vector_idx

【讨论】：

感谢@Morton 的精心回复！我感谢您付出的所有努力！我认为代码中存在一些错误，因为当我将返回的列元素相乘然后将点积与目标向量相乘时，分数与返回的分数不同。在第 16 行，您使用neg_matrix.sum(axis=1) 对每个 行 求和。这应该是axis=0 吗？实际上，返回的picked_vectors_idx 似乎是行索引而不是列索引，因为有时返回的值超出了列的范围。

啊，确实 - 我明白了。我假设向量的数量会大于它们的大小，但事实恰恰相反。我看看能不能改正。

好的，我相信你提到的实现是正确的。我忘了说我调换了你的表述。我发现使用行向量更容易一些。我只需要将矩阵形状更改为10,000x60,000 并将目标更改为60,000。我希望我早点意识到len(vec) > vec_num，因为由于丢弃了更多数据，所以计算速度要快得多。

该函数给出的分数仍然与我使用返回的索引计算它时计算的分数不同。我尝试使用返回的索引来索引MATRIX 中的行，但分数不同。在尝试他们对MATRIX 中的列进行索引时，我在计算点积时得到了形状不正确的错误。

我使用了您在代码顶部的示例输入数据、TARGET 和 MATRIX。我只是乘以元素返回的 MATRIX 的索引。然后我将这个结果与 TARGET 相乘。结果分数与函数返回的分数不同。

【参考方案2】：

也许是太天真了，但是我第一时间想到的就是选择离目标距离最短的5列：

import scipy
import numpy as np
from sklearn.metrics.pairwise import pairwise_distances

def sparse_prod_axis0(A):
    """Sparse equivalent of np.prod(arr, axis=0)
    From https://***.com/a/44321026/3381305
    """
    valid_mask = A.getnnz(axis=0) == A.shape[0]
    out = np.zeros(A.shape[1], dtype=A.dtype)
    out[valid_mask] = np.prod(A[:, valid_mask].A, axis=0)
    return np.matrix(out)

def get_strategy(M, target, n=5):
    """Guess n best vectors.
    """
    dists = np.squeeze(pairwise_distances(X=M, Y=target))
    idx = np.argsort(dists)[:n]
    return sparse_prod_axis0(M[idx])

# Example data.
M = scipy.sparse.rand(m=6000, n=1000, density=0.5, format='csr').astype('bool')
target = np.atleast_2d(np.random.choice([-1, 1], size=1000))

# Try it.
strategy = get_strategy(M, target, n=5)
result = strategy @ target.T

让我印象深刻的是，您可以再增加一个步骤，从 M–target 距离中取出前百分之几并检查它们的相互距离——但这可能会非常昂贵。

我没有检查这与详尽搜索相比如何。

【讨论】：

虽然这似乎是一种合乎逻辑的第一种方法，但向量的元素乘法完全改变了它们与目标的距离。例如，向量 c1 可能离目标很远，但是当与向量的正确组合相乘时，会产生最优策略 s。由于目标向量是稀疏的，我认为对于我的问题，最佳解决方案可能是密集向量与一个（最多两个稀疏）向量的混合。不幸的是，您的解决方案会产生五个稀疏向量（因为目标是稀疏的），当它们相乘时会导致向量 s 过于稀疏。

我想我应该能够利用这样的认识，即选择最多的 5 个向量应该是密集的，一个（最多两个）应该是稀疏的，但我还不太确定如何。

那么最大的和（按列）不应该起作用吗？我快速尝试了一下，但似乎没有帮助（将 dists 替换为 np.squeeze(np.asarray(np.sum(M, axis=0)))）。

strategy 给出了一个数组，其中所有内容都是“真”或“假”。我尝试使用 np.where(strategy)[1] 索引列号，并在您的示例数据上运行它时得到以下数组。 resulting_array = array([ 35, 52, 79, 355, 415, 417, 436, 495, 537, 549, 663, 686, 728, 737, 787, 814, 823, 832, 835, 925, 963])。难道我做错了什么？ resulting_array 中不应该是 5 个元素吗？谢谢！

似乎M[idx]（在第二个函数中返回）选择了五个 行 而不是列，因为它的形状是 5x1000，使用您的示例矩阵 M。我尝试将其更改为M[:,idx]，但随后它说idx 中的某些索引超出范围......似乎在get_strategy 函数中选择了行而不是列。 pairwise_distances 是否计算行之间的距离而不是列之间的距离可能是@kwinkunks？