<p dir="auto">Hello.<br>
I've got strange behavior for experiments. I'm working with matrix (for example <strong>b</strong>) that in result of multiplying <strong>b.T * b</strong> should be singular matrix and for inverse method should be arisen error like <em>numpy.linalg.linalg.LinAlgError: Singular matrix</em>. But result was high/low values.<br>
Code below:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; b = np.matrix([[1,1,0], [1,0,1], [1,1,0]])
&gt;&gt;&gt; b
matrix([[1, 1, 0],
        [1, 0, 1],
        [1, 1, 0]])
&gt;&gt;&gt; np.linalg.inv(b.T * b)
matrix([[ 4.50359963e+15, -4.50359963e+15, -4.50359963e+15],
        [-4.50359963e+15,  4.50359963e+15,  4.50359963e+15],
        [-4.50359963e+15,  4.50359963e+15,  4.50359963e+15]])"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; b = np.matrix([[1,1,0], [1,0,1], [1,1,0]])
&gt;&gt;&gt; b
matrix([[1, 1, 0],
        [1, 0, 1],
        [1, 1, 0]])
&gt;&gt;&gt; np.linalg.inv(b.T * b)
matrix([[ 4.50359963e+15, -4.50359963e+15, -4.50359963e+15],
        [-4.50359963e+15,  4.50359963e+15,  4.50359963e+15],
        [-4.50359963e+15,  4.50359963e+15,  4.50359963e+15]])
</code></pre></div>
<p dir="auto">How can be avoided this behavior?<br>
Tests on:<br>
win10, Python 3.5.4, numpy version '1.14.0'.<br>
ubuntu 16.04, Python 3.5.2, numpy version '1.13.3' and '1.14.0'.</p>
<p dir="auto">PS. I've checked via wolfram and R it's real singular matrix.</p>
<p dir="auto">I am trying to solve several independent systems of equations at the same time using <code class="notranslate">numpy.linalg.solve</code>, e.g., <code class="notranslate">a</code> has shape (N, M, M). The challenge I'm running into is how to deal with the case when, for some values of N, the last two dimensions comprise a singular matrix.</p>
<p dir="auto">Consider this example</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="import numpy as np
arr = [[[1, 0], [0, 1]], [[1, 0], [0, 1]], [[0, 0], [0, 0]]]
np.linalg.solve(arr, 1)"><pre class="notranslate"><span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>
<span class="pl-s1">arr</span> <span class="pl-c1">=</span> [[[<span class="pl-c1">1</span>, <span class="pl-c1">0</span>], [<span class="pl-c1">0</span>, <span class="pl-c1">1</span>]], [[<span class="pl-c1">1</span>, <span class="pl-c1">0</span>], [<span class="pl-c1">0</span>, <span class="pl-c1">1</span>]], [[<span class="pl-c1">0</span>, <span class="pl-c1">0</span>], [<span class="pl-c1">0</span>, <span class="pl-c1">0</span>]]]
<span class="pl-s1">np</span>.<span class="pl-s1">linalg</span>.<span class="pl-en">solve</span>(<span class="pl-s1">arr</span>, <span class="pl-c1">1</span>)</pre></div>
<p dir="auto">On numpy 1.9.2, this will raise a <code class="notranslate">LinAlgError</code>, but I still would like to know the solutions for all the full-rank matrices. (I would be okay with getting back NaNs for the singular cases.) I'm aware I can compute the singular values and use fancy indexing to slice the array, but in my algorithm I have to do several slicing/filtering steps and I would prefer not to lose the alignment of <code class="notranslate">a</code> with <code class="notranslate">b</code> -- my real-world <code class="notranslate">a</code> typically has 6-7 dimensions and all the indexing arrays make the code hard to follow.</p>
<p dir="auto">My ideal solution seems to involve masked arrays since I end up doing several other filtering steps, but it appears that most (all?) functions in the <code class="notranslate">linalg</code> family ignore array masks.</p>
<p dir="auto">Is there an alternative solution already available, or perhaps a suggestion on how I could contribute a solution that I could package as a PR?</p>