<p>这个问题有一个非常有效的(线性时间)解决方案,尽管它需要一些讨论。。。在</p>
<p><strong>Zeroth:澄清问题/LCP</strong></p>
<p>根据评论中的澄清,@FooBar说原始问题是元素方面的<code>min</code>;我们需要找到一个<code>z</code>(或<code>v</code>),这样</p>
<blockquote>
<p>either the left argument is zero and the right argument is non-negative or the left argument is non-negative and the right argument is zero</p>
</blockquote>
<p>正如@FooBar正确指出的,有效的重新参数化将导致LCP。然而,下面我展示了一个更简单、更有效的解决原始问题的方法(通过利用原始问题中的结构),而不需要LCP。为什么这应该更容易?注意,LCP在<code>z</code>(Bz+q)'z中有一个二次项,但是原始问题没有(只有线性项Bz+q和z)。下面我将利用这个事实。在</p>
<p><strong>首先:简化</strong></p>
<p>有一个重要但关键的细节在很大程度上简化了这个问题:</p>
<blockquote>
<ul>
<li>Create four matrices A11, A12, A21, A22 using scipy.sparse.diags </li>
<li>And
stack them together as A = scipy.sparse.bmat([[A11, A12], [A21,
A22]])</li>
</ul>
</blockquote>
<p>这有着巨大的影响。具体地说,这不是一个<strong>单一的</strong><strong>大的</strong>问题,而是一个<strong>非常小的</strong>问题(确切地说是2D)的<strong>大数问题。注意,这个<code>A</code>矩阵的块对角结构在所有后续操作中都保持不变。每一个后续的运算都是矩阵向量乘法或内积。这意味着这个程序实际上是在<code>z</code>(或<code>v</code>)变量对中的<em>可分离的</em>。在</p>
<p>具体来说,假设每个块<code>A11,...</code>的大小为<code>n</code>,大小为<code>n</code>。然后批判性地注意到<code>z_1</code>和{<cd12>}只在等式和术语中出现在</strong>中,而不会出现在其他地方。因此,该问题可分为<code>n</code>问题,每个问题都是二维的。因此,我将在以后解决2D问题,并且您可以在<code>n</code>上序列化或并行化,而不需要稀疏矩阵或大opt包。在</p>
<p><strong>第二:二维问题的几何结构</strong></p>
<p>假设我们有二维的元素问题,即:</p>
<pre><code>find z such that (elementwise) min( Bz + q , z ) = 0, or declare that no such `z` exists.
</code></pre>
<p>因为在我们的设置中<code>B</code>现在是2x2,这个问题几何对应于我们可以手动枚举的四个标量不等式(我将它们命名为a1、a2、z1、z2):</p>
^{pr2}$
<p>这代表一个(可能是空的)多面体,也就是二维空间中四个半空间的交集。在</p>
<p><strong>第三:解决二维问题</strong></p>
<p>(编辑:为了清楚起见,更新了这一点)</p>
<p>那么2D问题具体是什么呢?我们想要找到一个<code>z</code>,其中一个解决方案(虽然不是完全不同,但并不重要):</p>
<ol>
<li>a1>;=0,z1=0,a2>;=0,z2=0</li>
<li>a1=0,z1>;=0,a2=0,z2>;=0</li>
<li>a1>;=0,z1=0,a2=0,z2>;=0</li>
<li>a1=0,z1>;=0,a2>;=0,z2=0</li>
</ol>
<p>如果实现了其中之一,我们就有了一个解:z,其中z和Bz+q的元素最小值是0向量。如果这四个都不能实现,那么这个问题是不可行的,我们将宣布不存在这样的<code>z</code>。在</p>
<p>第一种情况发生在a1,a2的交点正或正中。只需选择解z=B^-1q,并检查它是否是元素非负的。如果是,则返回<code>B^-1q</code>(注意,即使B不是psd,我假设它是可逆的/满秩的)。所以:</p>
<pre><code>if B^-1q is elementwise nonnegative, return z = B^-1q.
</code></pre>
<p>第二种情况(与第一种情况不完全不同)发生在a1、a2的交点在任何地方但包含原点时。换句话说,当Bz+q>;0表示z=0时。当q为元素正时发生这种情况。所以:</p>
<pre><code>elif q is elementwise nonnegative, return z = 0.
</code></pre>
<p>第三种情况为z1=0,可以代入a2表示当z2=-q2/B22时a2=0。z2必须大于等于0,因此-q2/B22>;=0。a1也必须是>;=0,用z1和z2代替值,得到-B22/B12*q2+q1>;=0。所以:</p>
<pre><code>elif -q2/B22 >=0 and -B22/B12*q2 + q1 >=0, return z1= 0, z2 = -q2/B22.
</code></pre>
<p>第四步与第三步相同,但是交换1和2。所以:</p>
<pre><code>elif -q1/B11 >=0 and -B21/B11*q1 + q2 >=0, return z1 = -q1/B11, z2 =0.
</code></pre>
<p>最后第五种情况是当问题不可行时。当上述条件都不满足时会发生这种情况。所以:</p>
<pre><code>else return None
</code></pre>
<p><strong>最后:将各部分组合起来</strong></p>
<p>求解每一个二维问题都是一对简单/高效/平凡的线性解并进行比较。每个都将返回一对数字或<code>None</code>。然后对所有的<code>n</code>2D问题做同样的处理,并连接向量。如果有没有,这个问题是不可行的(全部没有)。否则,你有你的答案。在</p>