<p dir="auto">It's currently not possible to compute the full eigendecomposition with <code class="notranslate">eigsh</code>.</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; A = rand(5, 4)
&gt;&gt;&gt; AA = np.dot(A.T, A)
&gt;&gt;&gt; from scipy.sparse.linalg import eigsh
&gt;&gt;&gt; eigsh(AA, AA.shape[0])
[...]
ValueError: k must be between 1 and the order of the square input matrix."><pre class="notranslate"><span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-v">A</span> <span class="pl-c1">=</span> <span class="pl-en">rand</span>(<span class="pl-c1">5</span>, <span class="pl-c1">4</span>)
<span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-v">AA</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">dot</span>(<span class="pl-v">A</span>.<span class="pl-v">T</span>, <span class="pl-v">A</span>)
<span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-k">from</span> <span class="pl-s1">scipy</span>.<span class="pl-s1">sparse</span>.<span class="pl-s1">linalg</span> <span class="pl-k">import</span> <span class="pl-s1">eigsh</span>
<span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-en">eigsh</span>(<span class="pl-v">AA</span>, <span class="pl-v">AA</span>.<span class="pl-s1">shape</span>[<span class="pl-c1">0</span>])
[...]
<span class="pl-v">ValueError</span>: <span class="pl-s1">k</span> <span class="pl-s1">must</span> <span class="pl-s1">be</span> <span class="pl-s1">between</span> <span class="pl-c1">1</span> <span class="pl-c1">and</span> <span class="pl-s1">the</span> <span class="pl-s1">order</span> <span class="pl-s1">of</span> <span class="pl-s1">the</span> <span class="pl-s1">square</span> <span class="pl-s1">input</span> <span class="pl-s1">matrix</span>.</pre></div>
<p dir="auto"><a href="https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/arpack/arpack.py#L1506">https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/arpack/arpack.py#L1506</a></p>
<p dir="auto">I think <code class="notranslate">k &gt;= n</code> should be replaced by <code class="notranslate">k &gt; n</code>.</p>
<p dir="auto">In the context of <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1093032513" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15351" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/15351/hovercard" href="https://github.com/scipy/scipy/issues/15351">#15351</a> and reviewing <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1089902835" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15299" data-hovercard-type="pull_request" data-hovercard-url="/scipy/scipy/pull/15299/hovercard" href="https://github.com/scipy/scipy/pull/15299">#15299</a>, I noticed a bunch of inconsistencies that should be fixed IMO.</p>
<h2 dir="auto">TL;DR</h2>
<p dir="auto"><a href="https://github.com/scipy/scipy/issues/15600#issuecomment-1041258485" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/15600/hovercard">Courtesy</a> of <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/tirthasheshpatel/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/tirthasheshpatel">@tirthasheshpatel</a>:</p>
<ul dir="auto">
<li><code class="notranslate">factorial</code>, <code class="notranslate">factorial2</code>, and <code class="notranslate">factorialk</code> behave inconsistently for several different combinations of inputs. (for example, <code class="notranslate">factorial2</code> returning 0-dim arrays for scalars and, <code class="notranslate">factorial2</code> and <code class="notranslate">factorialk</code> aren't vectorized). See the first table for more details. This issue can be tacked by adding a new parameter like <code class="notranslate">legacy</code> which, by default, preserves the old behaviour but when false, handles the current inconsistent behavior correctly.</li>
<li>Support for float and complex inputs can easily be added using <a href="https://en.wikipedia.org/wiki/Double_factorial#Alternative_extension_of_the_multifactorial" rel="nofollow">this extension</a>. This enhancement adds a new boolean parameter <code class="notranslate">extend='zero'|'complex'</code> (and perhaps <code class="notranslate">'recurrence'</code> down the line) that specifies whether to use the extension formula or the default formula. See the second table for more details.</li>
</ul>
<h2 dir="auto">Design Space</h2>
<p dir="auto">There are several dimensions that come into play:</p>
<ul dir="auto">
<li>scalar vs. array inputs</li>
<li>for integers: approximation or exact computation (<code class="notranslate">exact=</code>)</li>
<li>extensions to negative integers</li>
<li>extensions to floats</li>
<li>extensions to complex numbers</li>
<li>extensions to poles (possible e.g. for <code class="notranslate">factorial2</code>)</li>
<li>for the array-case, failing soft (e.g. inserting NaNs) or hard (raising error)</li>
</ul>
<h2 dir="auto">Current Behaviour (as of 1.8.0)</h2>
<table role="table">
<thead>
<tr>
<th>Inputs</th>
<th><code class="notranslate">factorial(</code><br><code class="notranslate">exact=False)</code></th>
<th><code class="notranslate">factorial(</code><br><code class="notranslate">exact=True)</code></th>
<th><code class="notranslate">factorial2(</code><br><code class="notranslate">exact=False)</code></th>
<th><code class="notranslate">factorial2(</code><br><code class="notranslate">exact=True)</code></th>
<th><code class="notranslate">factorialk(</code><br><code class="notranslate">exact=False)</code></th>
<th><code class="notranslate">factorialk(</code><br><code class="notranslate">exact=True)</code></th>
</tr>
</thead>
<tbody>
<tr>
<td>Scalar positive integer</td>
<td>float</td>
<td>int</td>
<td>!! <em><strong>0-dim<br>float array</strong></em> !!</td>
<td>int</td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td>int</td>
</tr>
<tr>
<td>Scalar negative integer</td>
<td>0.0</td>
<td>0</td>
<td>!! <code class="notranslate">array(0.)</code> !!</td>
<td>0</td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td>0</td>
</tr>
<tr>
<td>Scalar positive float</td>
<td>float</td>
<td><code class="notranslate">ValueError</code> (but text misleading)</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Scalar negative float</td>
<td>0.0</td>
<td>0</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Scalar complex</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Scalar NaN</td>
<td>NaN</td>
<td>NaN</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Positive integer in array</td>
<td>float</td>
<td>int</td>
<td>float</td>
<td>!! <code class="notranslate">TypeError</code> !!</td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Negative integer in array</td>
<td>0.0</td>
<td>0</td>
<td>0.0</td>
<td>!! <code class="notranslate">TypeError</code> !!</td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Positive float in array</td>
<td>float</td>
<td><code class="notranslate">ValueError</code> (but text misleading)</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Negative float in array</td>
<td>0.0</td>
<td>0</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Complex in array</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>NaN in array</td>
<td>NaN</td>
<td>NaN</td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
<tr>
<td>Empty array</td>
<td><code class="notranslate">array([])</code></td>
<td><code class="notranslate">IndexError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">TypeError</code></td>
<td><code class="notranslate">NotImpl</code> <code class="notranslate">Error</code></td>
<td><code class="notranslate">TypeError</code></td>
</tr>
</tbody>
</table>
<p dir="auto">Note that:</p>
<ul dir="auto">
<li><code class="notranslate">factorial2(&lt;scalar&gt;, exact=False)</code> returns a zero-dimensional array instead of a scalar.</li>
<li>All the <code class="notranslate">TypeErrors</code> are IMO implementation oversights (and in case of "!!", against documentation) - i.e. we should provide better error messages.</li>
<li>There is a full <a href="https://en.wikipedia.org/wiki/Double_factorial#Alternative_extension_of_the_multifactorial" rel="nofollow">generalization</a> of <code class="notranslate">factorialk</code> which is basically just a fraction of two gamma functions, which would work for all complex n &amp; k (up to poles), as well as arrays.</li>
</ul>
<h2 dir="auto">Extension Issues</h2>
<p dir="auto">In addition to this hodgepodge of behaviour, some people would like to have the float/complex extensions of these functions (e..g. <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1092835847" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15349" data-hovercard-type="pull_request" data-hovercard-url="/scipy/scipy/pull/15349/hovercard" href="https://github.com/scipy/scipy/pull/15349">#15349</a>). I tend to agree, but that raises another problem:</p>
<h4 dir="auto">Complex extensions are incompatible with current inputs ≤0</h4>
<p dir="auto">First off, everything that's 0 currently on the negative axis gets a different value. But also note the "≤"; namely, 0!! is 1 in the regular definition, but ~0.79=sqrt(2/pi) in the <a href="https://en.wikipedia.org/wiki/Double_factorial#Complex_arguments" rel="nofollow">extension</a>.</p>
<h4 dir="auto">Extensions based on recurrence are incompatible with current status as well as complex extension</h4>
<p dir="auto">This came up for example in <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1089902835" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15299" data-hovercard-type="pull_request" data-hovercard-url="/scipy/scipy/pull/15299/hovercard" href="https://github.com/scipy/scipy/pull/15299">#15299</a>. It's <a href="https://en.wikipedia.org/wiki/Double_factorial#Negative_arguments" rel="nofollow">possible</a> (for example) to invert the recurrence relation n!! = n * (n-2)!!, and so extend the double-factorial to negative integers. The resulting value would be neither 0 (currently) nor the pole in the complex extension.</p>
<h2 dir="auto">Design Considerations</h2>
<p dir="auto">If we (or at least I) were designing from scratch:</p>
<ul dir="auto">
<li>Changing from the <code class="notranslate">exact=False</code> to <code class="notranslate">exact=True</code> should <em>never</em> fundamentally change the result (e.g. going from 0 to non-zero, or casting inputs to integers as seen in <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1093032513" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15351" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/15351/hovercard" href="https://github.com/scipy/scipy/issues/15351">#15351</a>).
<ul dir="auto">
<li>As such, the exact keyword only really matters to positive scalar integers; in all other cases, it could be retired, resp. raise, or simply be turned into a no-op (cf. <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1094046680" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15359" data-hovercard-type="pull_request" data-hovercard-url="/scipy/scipy/pull/15359/hovercard" href="https://github.com/scipy/scipy/pull/15359">#15359</a>)</li>
<li>Open question whether it's worth to properly implement <code class="notranslate">exact=True</code> for array-case (currently, only <code class="notranslate">factorial</code> does).</li>
</ul>
</li>
<li>Complex extensions  (together with non-positive floats) should be behind a keyword, roughly because:
<ul dir="auto">
<li>0-for-negative-values is a valid approach, and has been in scipy ~forever.</li>
<li>Since the values being produced are incompatible (for certain inputs), the complex extension should be explicitly opt-in.</li>
</ul>
</li>
<li>Using recurrence relations for extending into poles is an open question; not sure if worth it, but if so, probably with another <del>keyword</del> value for the <code class="notranslate">extend</code> keyword.</li>
<li>Use the most general approximation, except when alternatives with better accuracy / performance exist.</li>
<li>Probably need to think about dtype casting rules for empty arrays (keep dtype? always float? ...)</li>
</ul>
<h2 dir="auto">Possible Future Behaviour (note <code class="notranslate">extend=</code>, not <code class="notranslate">exact=</code>)</h2>
<table role="table">
<thead>
<tr>
<th>Inputs</th>
<th><code class="notranslate">factorial(</code><br><code class="notranslate">extend=</code><br><code class="notranslate">'zero')</code></th>
<th><code class="notranslate">factorial(</code><br><code class="notranslate">extend=</code><br><code class="notranslate">'complex')</code></th>
<th><code class="notranslate">factorial2(</code><br><code class="notranslate">extend=</code><br><code class="notranslate">'zero')</code></th>
<th><code class="notranslate">factorial2(</code><br><code class="notranslate">extend=</code><br><code class="notranslate">'complex')</code></th>
<th><code class="notranslate">factorialk(</code><br><code class="notranslate">extend=</code><br><code class="notranslate">'zero')</code></th>
<th><code class="notranslate">factorialk(</code><br><code class="notranslate">extend=</code><br><code class="notranslate">'complex')</code></th>
</tr>
</thead>
<tbody>
<tr>
<td>Scalar positive integer w/ <code class="notranslate">exact=True</code></td>
<td>n!</td>
<td>n!</td>
<td>n!!</td>
<td>n!!</td>
<td>n!(k)</td>
<td>n!(k)</td>
</tr>
<tr>
<td>Scalar positive integer</td>
<td>Γ(n+1)</td>
<td>Γ(n+1)</td>
<td>(1)</td>
<td>(1)</td>
<td>(2)</td>
<td>(2)</td>
</tr>
<tr>
<td>Scalar negative integer / poles</td>
<td>0.0</td>
<td>NaN</td>
<td>0.0</td>
<td>NaN</td>
<td>0.0</td>
<td>NaN</td>
</tr>
<tr>
<td>Scalar positive float</td>
<td>Γ(n+1)</td>
<td>Γ(n+1)</td>
<td>(1)</td>
<td>(1)</td>
<td>(2)</td>
<td>(2)</td>
</tr>
<tr>
<td>Scalar negative float</td>
<td>0.0</td>
<td>Γ(n+1)</td>
<td>0.0</td>
<td>(1)</td>
<td>0.0</td>
<td>(2)</td>
</tr>
<tr>
<td>Scalar complex</td>
<td><code class="notranslate">ValueError</code></td>
<td>Γ(n+1)</td>
<td><code class="notranslate">ValueError</code></td>
<td>(1)</td>
<td><code class="notranslate">ValueError</code></td>
<td>(2)</td>
</tr>
<tr>
<td>Scalar NaN</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
</tr>
<tr>
<td>Positive integer in array</td>
<td>Γ(n+1)</td>
<td>Γ(n+1)</td>
<td>(1)</td>
<td>(1)</td>
<td>(2)</td>
<td>(2)</td>
</tr>
<tr>
<td>Negative integer in array</td>
<td>0.0</td>
<td>NaN</td>
<td>0.0</td>
<td>NaN</td>
<td>0.0</td>
<td>NaN</td>
</tr>
<tr>
<td>Positive float in array</td>
<td>Γ(n+1)</td>
<td>Γ(n+1)</td>
<td>(1)</td>
<td>(1)</td>
<td>(2)</td>
<td>(2)</td>
</tr>
<tr>
<td>Negative float in array</td>
<td>0.0</td>
<td>Γ(n+1)</td>
<td>0.0</td>
<td>(1)</td>
<td>0.0</td>
<td>(2)</td>
</tr>
<tr>
<td>Complex in array</td>
<td><code class="notranslate">ValueError</code></td>
<td>Γ(n+1)</td>
<td><code class="notranslate">ValueError</code></td>
<td>(1)</td>
<td><code class="notranslate">ValueError</code></td>
<td>(2)</td>
</tr>
<tr>
<td>NaN in array</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
<td>NaN</td>
</tr>
<tr>
<td>Empty array</td>
<td><code class="notranslate">array([])</code></td>
<td><code class="notranslate">array([])</code></td>
<td><code class="notranslate">array([])</code></td>
<td><code class="notranslate">array([])</code></td>
<td><code class="notranslate">array([])</code></td>
<td><code class="notranslate">array([])</code></td>
</tr>
</tbody>
</table>
<p dir="auto">The formulas in question are:</p>
<ul dir="auto">
<li>(1): z!! = 2^((z-1)/2) Γ(z/2 + 1) / Γ(1/2 + 1)</li>
<li>(2): z!(α) =  α^((z-1)/α) Γ(z/α + 1) / Γ(1/α + 1)</li>
</ul>
<p dir="auto">I'd very much welcome feedback on this - there are some trade-offs to be made (e.g. whether negative floats should also raises without <code class="notranslate">extend='complex'</code> or whether the errors in the array-case should just map to NaN; also there are different approximations for n!! in use <a href="https://github.com/scipy/scipy/blob/v1.8.0/scipy/special/_basic.py#L2411-L2412">currently</a>), but I think the above would be reasonable. In particular, the scalar and array case behave the same way.</p>
<h2 dir="auto">Implementation Plan</h2>
<p dir="auto">I would propose to introduce another keyword (sorry...) <code class="notranslate">legacy</code>, defaulting to <code class="notranslate">True</code>, but implement <code class="notranslate">legacy=False</code> with all the corrected behaviour. Where there's any planned change in behaviour, we raises a deprecation warning recommending to switch to <code class="notranslate">legacy=False</code>. After 2 releases, we switch the default from <code class="notranslate">legacy=True</code> to <code class="notranslate">legacy=False</code>, after two more releases we raise on <code class="notranslate">legacy=True</code>, and after two more we remove the keyword.</p>
<p dir="auto">Looking more closely at relevant changes in (currently non-raising) behaviour, it really only affects the zero-dimensional arrays for <code class="notranslate">factorial2(scalar_val, exact=False)</code>. In particular, if we let <code class="notranslate">exact=True</code> be a no-op for non-integers (the less "noisy" options from the point of UX), we would keep essentially all current non-raising behaviour intact.</p>
<p dir="auto">CC <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/steppi/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/steppi">@steppi</a> <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/tirthasheshpatel/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/tirthasheshpatel">@tirthasheshpatel</a> <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/j-bowhay/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/j-bowhay">@j-bowhay</a> <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/tupui/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/tupui">@tupui</a></p>
<h3 dir="auto">Edit History:</h3>
<ul dir="auto">
<li>v4 (2022-03-10): Make <code class="notranslate">exact=True</code> a no-op for non-integers (less noisy), following <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1094046680" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15359" data-hovercard-type="pull_request" data-hovercard-url="/scipy/scipy/pull/15359/hovercard" href="https://github.com/scipy/scipy/pull/15359">#15359</a>. Extend deprecation period of <code class="notranslate">legacy</code>-keyword along the same lines as proposed in <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="1164529699" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/15722" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/15722/hovercard" href="https://github.com/scipy/scipy/issues/15722">#15722</a>.</li>
<li>v3 (2022-02-16; up to rev 20): Added (current &amp; proposed) behaviour for NaN and empty arrays; good catch <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/tirthasheshpatel/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/tirthasheshpatel">@tirthasheshpatel</a>!</li>
<li>v2 (2022-02-16; up to rev 15): Added TL;DR, clarification on impact of <code class="notranslate">legacy</code>, changed <code class="notranslate">extend</code> from boolean argument to take an enum.</li>
<li>v1 (2022-02-16; up to rev 10): Original</li>
</ul>