最小差约束稀疏最小二乘问题（从 python 调用）答案

【问题标题】：minimum-difference constrained sparse least squares problem (called from python)最小差约束稀疏最小二乘问题（从 python 调用）
【发布时间】：2022-11-13 12:18:08
【问题描述】：

我正在努力寻找适合的快速算法。

我只想最小化：

范数2(x-s)

英石 G.x <= h

x >= 0

总和（x）= R

G 是稀疏的并且仅包含 1（显然还有零）。

在迭代算法的情况下，最好将临时解决方案展示给用户。

上下文是 s 是当前结果的向量，用户说“这几个条目的总和（由 G 中的一行中的几个 1.0 表示的条目）应该小于这个值（在h). 所以我们必须以最小二乘最优方式从用户指定的条目（由 G 中的 1.0 条目表示）中删除数量，但是由于我们对总数 (R) 有一个全局约束，因此删除的值需要是在其他条目中以最小二乘最佳方式分配。条目不能为负数。

我正在查看的所有算法都是很多更一般，因此也更复杂。而且，它们似乎很慢。我不认为这是一个复杂的问题，尽管平等和不平等约束的混合似乎总是让事情变得更加复杂。

这必须从 Python 调用，所以我正在查看 qpsolvers 和 scipy.optimize 等 Python 库。但我想 Java 或 C++ 库可以从 Python 中使用和调用，这可能很好，因为多线程在 Java 和 C++ 中更好。

关于使用什么库/包/方法来最好地解决这个问题有什么想法吗？

问题的大小在 s 中大约有 150,000 行，在 G 中只有几十行。

谢谢！

【问题讨论】：

标签： least-squares scipy-optimize

【解决方案1】：

您的问题是线性最小二乘法：

minimize_x   norm2(x-s)
such that    G x <= h
             x >= 0
             1^T x = R

它符合 qpsolvers 中 solve_ls 函数的要求。

鉴于您指定的内容，这是我想象您的问题矩阵想要的一个实例。由于它是稀疏的，我们应该确保使用 SciPy CSC 矩阵（以及用于向量的常规 NumPy 数组）：

import numpy as np
import scipy.sparse as spa

n = 150_000

# minimize || x - s ||^2
R = spa.eye(n, format="csc")
s = np.array(range(n), dtype=float)

# such that G * x <= h
G = spa.diags(
    diagonals=[
        [1.0 if i % 2 == 0 else 0.0 for i in range(n)],
        [1.0 if i % 3 == 0 else 0.0 for i in range(n - 1)],
        [1.0 if i % 5 == 0 else 0.0 for i in range(n - 1)],
    ],
    offsets=[0, 1, -1],
)
a_dozen_rows = np.linspace(0, n - 1, 12, dtype=int)
G = G[a_dozen_rows]
h = np.ones(12)

# such that sum(x) == 42
A = spa.csc_matrix(np.ones((1, n)))
b = np.array([42.0]).reshape((1,))

# such that x >= 0
lb = np.zeros(n)

接下来，我们可以解决这个问题：

from qpsolvers import solve_ls

x = solve_ls(R, s, G, h, A, b, lb, solver="cvxopt", verbose=True)

在这里，我选择了 CVXOPT，但您还可以安装其他开源求解器，例如 ProxQP、OSQP 或 SCS。您可以通过以下方式安装一组开源求解器：pip install qpsolvers[open_source_solvers]。安装一些求解器后，您可以通过以下方式列出稀疏矩阵的求解器：

print(qpsolvers.sparse_solvers)

最后，这里有一些代码来检查求解器返回的解决方案是否满足我们的约束：

tol = 1e-8  # tolerance for checks
print(f"- Objective: {np.linalg.norm(x - s):.1f}")
print(f"- G * x <= h: {(G.dot(x) <= h + tol).all()}")
print(f"- x >= 0: {(x + tol >= 0.0).all()}")
print(f"- sum(x) = {x.sum():.1f}")

希望这会有所帮助。快乐解决！

【讨论】：