Repository: scipy/scipy
Pull request: BUG: linalg: support empty arrays
#### Reference issue
<!--Example: Closes gh-WXYZ.-->

supersedes and closes gh-18673, supersedes and closes gh-18141

closes https://github.com/scipy/scipy/issues/17658
closes https://github.com/scipy/scipy/issues/14389
closes https://github.com/scipy/scipy/issues/14244
closes https://github.com/scipy/scipy/issues/9459
closes https://github.com/scipy/scipy/issues/8510
closes https://github.com/scipy/scipy/issues/8056

#### What does this implement/fix?
<!--Please explain your changes.-->

Includes two previous PRs, squashed to simplify the rebasing.
An annoying thing is making sure the dtype is right: `solve(int, float32) -> float64`, `solve(float32, float32) -> float32` etc.  
We absolutely need to make sure that empty and non-empty code paths agree.

Functions (`dir(linalg)`, slightly edited):

- [ ] bandwidth
- [ ] block_diag
- [ ] cdf2rdf
- [x] cho_factor
- [x] cho_solve
- [x] cho_solve_banded
- [x] cholesky
- [x] cholesky_banded
- [ ] circulant
- [ ] clarkson_woodruff_transform
- [ ] companion
- [ ] convolution_matrix
- [ ] coshm
- [ ] cosm
- [ ] cossin
- [x] det
- [ ] dft
- [ ] diagsvd
- [x] eig
- [x] eig_banded
- [x] eigh
- [x] eigh_tridiagonal
- [x] eigvals
- [x] eigvals_banded
- [x] eigvalsh
- [x] eigvalsh_tridiagonal
- [x] expm
- [ ] expm_cond
- [ ] expm_frechet
- [ ] fiedler
- [ ] fiedler_companion
- [ ] fractional_matrix_power
- [ ] funm
- [ ] hadamard
- [ ] hankel
- [ ] helmert
- [x] hessenberg
- [ ] hilbert
- [x] inv
- [ ] invhilbert
- [ ] invpascal
- [ ] khatri_rao
- [ ] kron
- [ ] ldl
- [ ] leslie
- [x] logm
- [x] lstsq
- [x] lu
- [x] lu_factor
- [x] lu_solve
- [ ] matmul_toeplitz
- [x] matrix_balance
- [ ] norm
- [x] null_space
- [ ] ordqz
- [x] orth
- [ ] orthogonal_procrustes
- [ ] pascal
- [x] pinv
- [x] pinvh
- [ ] polar
- [x] qr
- [ ] qr_delete
- [ ] qr_insert
- [x] qr_multiply
- [ ] qr_update
- [ ] qz
- [x] rq
- [ ] rsf2csf
- [x] schur
- [ ] signm
- [ ] sinhm
- [ ] sinm
- [x] solve
- [x] solve_banded
- [x] solve_circulant
- [ ] solve_continuous_are
- [ ] solve_continuous_lyapunov
- [ ] solve_discrete_are
- [ ] solve_discrete_lyapunov
- [ ] solve_lyapunov
- [ ] solve_sylvester
- [x] solve_toeplitz
- [x] solve_triangular
- [x] solveh_banded
- [x] sqrtm
- [ ] subspace_angles
- [x] svd
- [x] svdvals
- [ ] tanhm
- [ ] tanm
- [ ] toeplitz


#### Additional information
<!--Any additional information you think is important.-->

Would be great to have a review of the main approach in the last commits --- `solve`, `solve_triangular` and `inv` before I manually go over the rest of linalg functions

cc @JozsefKutas @boatwrong 


Problem statement:
cho_factor and cho_solve don't support (0,0)-shape matrices
`scipy.linalg.cho_factor` doesn't support arrays with shape `(0, 0)`. I think it could be supported and this would avoid writing if-else conditions outside this function.

### Reproducing code example:
```python
import scipy
scipy.linalg.cho_factor(np.ones((0,0)))
```

### Error message:
```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-35b2e5bb3eee> in <module>()
----> 1 scipy.linalg.cho_factor(np.ones((0,0)))

/.../lib/python3.5/site-packages/scipy/linalg/decomp_cholesky.py in cho_factor(a, lower, overwrite_a, check_finite)
    130     """
    131     c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=False,
--> 132                             check_finite=check_finite)
    133     return c, lower
    134 

/.../lib/python3.5/site-packages/scipy/linalg/decomp_cholesky.py in _cholesky(a, lower, overwrite_a, clean, check_finite)
     26     overwrite_a = overwrite_a or _datacopied(a1, a)
     27     potrf, = get_lapack_funcs(('potrf',), (a1,))
---> 28     c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean)
     29     if info > 0:
     30         raise LinAlgError("%d-th leading minor not positive definite" % info)

ValueError: On entry to DPOTRF parameter number 4 had an illegal value
```

### Scipy/Numpy/Python version information:
```
0.18.1 1.11.1 sys.version_info(major=3, minor=5, micro=3, releaselevel='final', serial=0)
```



ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)
I tried to perform matrix inversion using 'scipy.linalg.inv', then the error below emerged : 
  File "D:\Anaconda2\lib\site-packages\scipy\linalg\basic.py", line 806, in inv
    lu, piv, info = getrf(a1, overwrite_a=overwrite_a)

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)
Does anyone know how this happened? Thanks!

### Reproducing code example:
numpy.linalg.inv( B.T.dot( B ) + ( 2 * gamma / penalty_coeff ) * numpy.identity( B.shape[ 1 ], dtype = numpy.float64 ) )
### Error message:
ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)

### Scipy/Numpy/Python version information:
python version : 2.7.12 from Anaconda 4.1.1 ( 64-bit )
numpy version : 1.13.1
scipy.linalg version : 0.4.9


BUG: linalg: lu and decompositions don't support (0, 1) or (0, 0) shape input
<!-- 
Thank you for taking the time to report a SciPy issue.
Please describe the issue in detail, and for bug reports
fill in the fields below. You can delete the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
-->

LU and QR decompositions return correct output matrices if input matrix shape is (1, 0), but fail for (0, 0) or (0, 1) shaped inputs.

### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be marked as a code block automatically 
-->


```
import numpy as np
from scipy import linalg

linalg.qr(np.random.randn(1, 0))  # Returns expected empty output matrices
linalg.qr(np.random.randn(0, 1))
```

### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
(array([[1.]]), array([], shape=(1, 0), dtype=float64))

ValueErrorTraceback (most recent call last)
<ipython-input-320-acec9ee93167> in <module>()
      3 
      4 print(linalg.qr(np.random.randn(1, 0)))  # Returns expected empty output matrices
----> 5 print(linalg.qr(np.random.randn(0, 1)))

/usr/local/python/python-2.7/std/lib/python2.7/site-packages/scipy/linalg/decomp_qr.py in qr(a, overwrite_a, lwork, mode, pivoting, check_finite)
    139         geqrf, = get_lapack_funcs(('geqrf',), (a1,))
    140         qr, tau = safecall(geqrf, "geqrf", a1, lwork=lwork,
--> 141                            overwrite_a=overwrite_a)
    142 
    143     if mode not in ['economic', 'raw'] or M < N:

/usr/local/python/python-2.7/std/lib/python2.7/site-packages/scipy/linalg/decomp_qr.py in safecall(f, name, *args, **kwargs)
     17     if lwork in (None, -1):
     18         kwargs['lwork'] = -1
---> 19         ret = f(*args, **kwargs)
     20         kwargs['lwork'] = ret[-2][0].real.astype(numpy.int)
     21     ret = f(*args, **kwargs)

ValueError: On entry to Dgeqrf parameter number 4 had an illegal value
```

### Scipy/Numpy/Python version information:
```
'1.1.0', '1.14.3', sys.version_info(major=2, minor=7, micro=15, releaselevel='final', serial=0)
```


ValueError: On entry to DGESDD parameter number 10 had an illegal value
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

My issue is about ... the library scipy, I used the classifier : LDA to train on a small dataset but I keep getting the error below : 
ValueError: On entry to DGESDD parameter number 10 had an illegal value

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```

Sample code to reproduce the problem

# scikit-learn==0.22.1
# pandas == 1.0.3
# numpy == 1.18.4
# scipy == 1.4.1
# torch == 1.5.0a0
# python == 3.7.6
# matplotlib==3.2.2
# transformers==4.4.2


model8 = XGBClassifier(eval_metric="logloss")
model7 = LinearDiscriminantAnalysis()
models_list = {'Model_XGB': model8, 'Model_LDA': model7}

# cross_validation

cv_splitter = KFold(n_splits=10, shuffle=True, random_state=42)
features, s = get_flaubert_layer(features, lge_model)
#print(features)
for model_name, model in models_list.items():
print("Model name: {}".format(model_name))
print("------------cross val predict used----------------", "\n")
y_pred = cross_val_predict(model, features, ylabels, cv=cv_splitter, verbose=1)
accuracy_score_predict = accuracy_score(ylabels, y_pred)
print("accuracy with cross_val_predict :", accuracy_score_predict)


```

#### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  Loading pytorch-gpu/py3/1.5.0
  Loading requirement: gcc/8.3.1 cuda/10.1.2 nccl/2.5.6-2-cuda
    cudnn/10.1-v7.5.1.10 intel-compilers/19.0.4 openmpi/4.0.1-cuda
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:    2.7s finished
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:    3.6s finished
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:    1.7s finished
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:  1.3min finished
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed: 15.6min finished
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:   47.8s finished
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
Traceback (most recent call last):
  File "training_for_multiclasse_classification.py", line 62, in <module>
    main()
  File "training_for_multiclasse_classification.py", line 44, in main
    dic_acc, dic_report, dic_cm, s = cross_validation(data_train, data_label_train, models_list, name, language_model_dir)
  File "../traitements/functions_for_training.py", line 58, in cross_validation
    y_pred = cross_val_predict(model, features, ylabels, cv=cv_splitter, verbose=1)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/sklearn/model_selection/_validation.py", line 755, in cross_val_predict
    for train, test in cv.split(X, y, groups))
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/joblib/parallel.py", line 1029, in __call__
    if self.dispatch_one_batch(iterator):
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/joblib/parallel.py", line 847, in dispatch_one_batch
    self._dispatch(tasks)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/joblib/parallel.py", line 765, in _dispatch
    job = self._backend.apply_async(batch, callback=cb)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/joblib/_parallel_backends.py", line 206, in apply_async
    result = ImmediateResult(func)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/joblib/_parallel_backends.py", line 570, in __init__
    self.results = batch()
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/joblib/parallel.py", line 253, in __call__
    for func, args, kwargs in self.items]
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/joblib/parallel.py", line 253, in <listcomp>
    for func, args, kwargs in self.items]
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/sklearn/model_selection/_validation.py", line 841, in _fit_and_predict
    estimator.fit(X_train, y_train, **fit_params)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/sklearn/discriminant_analysis.py", line 477, in fit
    self._solve_svd(X, y)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/sklearn/discriminant_analysis.py", line 396, in _solve_svd
    _, S, V = linalg.svd(X, full_matrices=0)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/scipy/linalg/decomp_svd.py", line 125, in svd
    compute_uv=compute_uv, full_matrices=full_matrices)
  File "/gpfslocalsup/pub/anaconda-py3/2020.02/envs/pytorch-gpu-1.5.0/lib/python3.7/site-packages/scipy/linalg/lapack.py", line 806, in _compute_lwork
    ret = routine(*args, **kwargs)
ValueError: On entry to DGESDD parameter number 10 had an illegal value
srun: error: r6i7n6: task 0: Exited with exit code 1
srun: Terminating job step 221450.0


# X equal 


[[ 0.  0. -0. ... -0. -0. -0.]
 [ 0. -0. -0. ... -0. -0. -0.]
 [ 0.  0.  0. ... -0. -0.  0.]
 ...
 [ 0. -0. -0. ... -0. -0. -0.]
 [ 0. -0. -0. ... -0. -0. -0.]
 [ 0.  0. -0. ... -0. -0. -0.]]

```

```

#### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
--> 
1.4.1 1.18.4 sys.version_info(major=3, minor=7, micro=6, releaselevel='final', serial=0)




`linalg.inv` fails for arrays of shape (0, 0)
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

My issue is about the behavior of `linalg.inv` for arrays of shape (0, 0).
I would like to move from `numpy.linalg` to `scipy.linalg`, but with scipy inverting this kind of matrices is not possible.
Normally my functions runs with (M, M) matrices with M > 0 but I would like to handle the edge case as well where M = 0.
I can work around this issue by checking the shape beforehand, but it seams to me like a bug with scipy's `linalg.inv` that it can't handle this case.

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```
import numpy
import scipy
z = numpy.zeros((0, 0))
numpy.linalg.inv(z) # --> array([], shape=(0, 0), dtype=float64)
scipy.linalg.inv(z) # --> ValueError
```

#### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/test.py in <module>
      3 z = numpy.zeros((0, 0))
      4 numpy.linalg.inv(z) # --> array([], shape=(0, 0), dtype=float64)
----> 5 scipy.linalg.inv(z) # --> ValueError

~/anaconda3/lib/python3.8/site-packages/scipy/linalg/basic.py in inv(a, overwrite_a, check_finite)
    963         raise LinAlgError("singular matrix")
    964     if info < 0:
--> 965         raise ValueError('illegal value in %d-th argument of internal '
    966                          'getrf|getri' % -info)
    967     return inv_a

ValueError: illegal value in 4-th argument of internal getrf|getri
```

#### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
-->
`1.6.2 1.20.2 sys.version_info(major=3, minor=8, micro=5, releaselevel='final', serial=0)`

Inconsistent support for empty matrices in linalg
Some functions in `scipy.linalg` handle empty matrix arguments and some don't - there doesn't seem to be an obvious pattern. On the other hand, `numpy.linalg` seems to support empty matrices in all functions (though I didn't test exhaustively).

For example, `solve` and `cholesky` both support empty matrices by explicitly handling empty matrix arguments before calling LAPACK:

```python
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([]).reshape((0, 0))
>>> linalg.solve(a, a)
array([], shape=(0, 0), dtype=float64)
>>>  linalg.cholesky(a)
array([], shape=(0, 0), dtype=float64)
```

while `inv` and most other functions don't:

```python
>>> linalg.inv(a)
```

```
Traceback (most recent call last):
  File "...\.venv\Lib\site-packages\IPython\core\interactiveshell.py", line 3433, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-12-8c6d9a341a08>", line 1, in <module>
    linalg.inv(a)
  File "...\.venv\Lib\site-packages\scipy\linalg\_basic.py", line 958, in inv
    raise ValueError('illegal value in %d-th argument of internal '
ValueError: illegal value in 4-th argument of internal getrf|getri
```

Empty matrix support can be convenient when matrix sizes are determined dynamically, for example, according to statistical significance tests.

Hints from issue discussion:
And similarly for `cho_solve`.
Seems to work in 1.0.0rc2

... except cho_solve apparently doesn't

This is a bit unfortunate that these are exposed to the public API. `cho_solve` only works with the input produced by `cho_factor` and for this special case it is not guarded.

In the ideal world there would be only a cholesky solver which combines both these functions. Should I attempt to do that or is it not worth deprecating them? Because `cholesky` already does the factorization if only factorization needed.
Same story with `lu_factor` and `lu_solve`. This is already in the `roadmap.rst`
About deprecating cho_solve:

Is there an alternative now for repeated cho_solve using the same "A" array?

I haven't looked at this in years, but the usecase for this in linalg and sparse.linalg was to avoid repeated computation of the cho or lu factorization.

@josef-pkt `cho_solve` is literally thin wrapper around `?potrs` with size, `check_finite` and `_asarray` checks. So you can basically change the function name. That's why I'm a bit skeptical about their exposition to the public API. 

Of course the counter argument is forcing users to use a low-level LAPACK routine which is also not ideal.  However, I think we could do better with the API. 
We use cho_solve in a few places in statsmodels, only one is recent, AFAIK,

I like thin wrappers if they do the tedious parts and we don't have to remember the low level interface nor the LAPACK names.
I should clarify, imagine there were only cholesky and solve_cholesky in the public API and solve_cholesky had a signature like 

    def solve_cholesky(a, b, factored_a, lower.....)
        """
        ....
        factored_a : (n, n) array_like
            The cholesky-factored coefficient matrix to bypass the factorization for quick reuse.
        ....
        """

Well of course wording can be improved but all are in one place and the distinction is clear and avoids duplication of cholesky and cho_factor. Also fits in well with the other function names in the linalg module. 
> def solve_cholesky(a, b, factored_a, lower.....)

but you don't need `a` anymore, so the signature becomes more messy than just the plain chol_solve
(unless a would autodetect that it is actually factored_a
That should not have happened but still it might happen if you supplied an empty matrix to inversion function. What is the resulting array? 
Sorry, I've made an mistake : the code should be 'scipy.linalg.inv(...' ; I guess it's more probable that the error has something to do with 'scipy.linalg.inv', since no error was raised when I re-ran the code with 'numpy.linalg.inv' in its place.
Ok what is the resulting array and your version of SciPy? 
I thought there was no array returned since the error raised should have prevented the execution. I noticed that 'numpy.linalg.inv' will return an emty array if it's fed with an empty 2-D array, instead of the 'ValueError' here.
Here's a simple example that produces the error by passing an array with shape (0,0) to `scipy.linalg.inv`:

```
In [2]: import numpy as np

In [3]: import scipy.linalg as la

In [4]: a = np.array([]).reshape(0, 0)

In [5]: a
Out[5]: array([], shape=(0, 0), dtype=float64)

In [6]: la.inv(a)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-543afa477a92> in <module>()
----> 1 la.inv(a)

~/miniconda3scipy/lib/python3.5/site-packages/scipy-1.1.0.dev0+a564e8f-py3.5-macosx-10.6-x86_64.egg/scipy/linalg/basic.py in inv(a, overwrite_a, check_finite)
    960                                                   'getri_lwork'),
    961                                                  (a1,))
--> 962     lu, piv, info = getrf(a1, overwrite_a=overwrite_a)
    963     if info == 0:
    964         lwork = _compute_lwork(getri_lwork, a1.shape[0])

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)
```
@WarrenWeckesser You are right. The reseponse of 'scipy.linalg.inv' and that of 'numpy.linalg.inv' are different when both are fed with the same empty 2-D array.
@HawkRong What is your matrix that leads to the error? Do you also have an empty array? 
Hi @keloemma Could you please identify exactly the value that hits `svd`  on this line in the stacktrace? 

```
_, S, V = linalg.svd(X, full_matrices=0)
```

Otherwise we don't know what happens internally inside the scikit-learn part and probably should be better suited for sklearn folks to have a look. Seems like at some point this `X`  becomes an empty array. 


@ilayn  When printing X , I get an array of 0 elements. (see below)
Thank you for raising this. It is indeed inconsistent with what NumPy does. @ilayn what do you think? Should we improve the validation here?
Argh. Yes we need to accomodate for these indeed. 

But as a word of obligatory caution, prefer `solve` to `inv` if you don't need the inverse itself. 
@Kaniee would you like to make a PR?
Would be my first one here, but I can give it at try.
I would add something like this:
```
if a1.shape[0] == 0:
    return np.identity(0)
```
You can just return `a1.copy()` after the squareness check. 
There are several open issues about zero-sized matrices in other linalg functions, would be nice to handle these uniformly across scipy.linalg.

IIRC there also was an issue with a link to a Nick Higham paper about zero-sized matrices, but I cannot find it anymore. IIRC it was Alex Griffing advocating for it. 
@ev-br  #4908 ? 
Thanks @ilayn, clicking though... found it ---It's de Boor, not Higham, it turns out:
https://github.com/scipy/scipy/pull/4672#issuecomment-124718692 
Agreed. This is something we should consider breaking backwards compatibility to attain a consistent implementation. Edge cases had been our soft spot in linalg for such functions that grew a bit heterogeneously. 
I'm interested in looking into this issue. From what I can tell the cholesky decomp function returns an empty array of the same size. This seems like it would also serve as a valid return for the inverse function right? 

Also, I'll look into other functions that could benefit from empty matrix support that do not already implement it. 
> I'm interested in looking into this issue. From what I can tell the cholesky decomp function returns an empty array of the same size. This seems like it would also serve as a valid return for the inverse function right?

Yup, it should just be an array of the same (zero size) dimensions as the input.

I think either `solve` or `cholesky` could be used for reference. For functions like `inv`, the result should be consistent with the result from the corresponding function in numpy.

Hi @boatwrong, thanks for working on this. I should have some free time over the next couple of weeks, so if you want I'd be happy to continue from where you left off.
@JosefKutas Apologies for my late response, thank you for reaching out and looking into this! Please feel free to continue on.

What you’ve explained/written looks great and I appreciate you taking it further than I could.
________________________________
From: Jozsef Kutas ***@***.***>
Sent: Wednesday, June 7, 2023 18:23
To: scipy/scipy ***@***.***>
Cc: boatwrong ***@***.***>; Mention ***@***.***>
Subject: Re: [scipy/scipy] Inconsistent support for empty matrices in linalg (Issue #17658)


Hi @boatwrong<https://github.com/boatwrong>, thanks for working on this. I should have some free time over the next couple of weeks, so if you want I'd be happy to continue from where you left off.

—
Reply to this email directly, view it on GitHub<https://github.com/scipy/scipy/issues/17658#issuecomment-1581648098>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AUUREPEXPWJU422HRC7XIP3XKEENLANCNFSM6AAAAAATJ6YOQU>.
You are receiving this because you were mentioned.Message ID: ***@***.***>

Commit message: MAINT: fix issue with undefined variable in linalg.svd

Introduced because of a cross-merge conflict with another PR.

Patch to review:
diff --git a/scipy/linalg/_basic.py b/scipy/linalg/_basic.py
index c781236f2573..1fdea328693c 100644
--- a/scipy/linalg/_basic.py
+++ b/scipy/linalg/_basic.py
@@ -163,7 +163,8 @@ def solve(a, b, lower=False, overwrite_a=False,
 
     # accommodate empty arrays
     if b1.size == 0:
-        return np.asfortranarray(b1.copy())
+        dt = solve(np.eye(2, dtype=a1.dtype), np.ones(2, dtype=b1.dtype)).dtype
+        return np.empty_like(b1, dtype=dt)
 
     # regularize 1-D b arrays to 2D
     if b1.ndim == 1:
@@ -332,10 +333,20 @@ def solve_triangular(a, b, trans=0, lower=False, unit_diagonal=False,
 
     a1 = _asarray_validated(a, check_finite=check_finite)
     b1 = _asarray_validated(b, check_finite=check_finite)
+
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
         raise ValueError('expected square matrix')
+
     if a1.shape[0] != b1.shape[0]:
         raise ValueError(f'shapes of a {a1.shape} and b {b1.shape} are incompatible')
+
+    # accommodate empty arrays
+    if b1.size == 0:
+        dt_nonempty = solve_triangular(
+            np.eye(2, dtype=a1.dtype), np.ones(2, dtype=b1.dtype)
+        ).dtype
+        return np.empty_like(b1, dtype=dt_nonempty)
+
     overwrite_b = overwrite_b or _datacopied(b1, b)
 
     trans = {'N': 0, 'T': 1, 'C': 2}.get(trans, trans)
@@ -429,15 +440,22 @@ def solve_banded(l_and_u, ab, b, overwrite_ab=False, overwrite_b=False,
 
     a1 = _asarray_validated(ab, check_finite=check_finite, as_inexact=True)
     b1 = _asarray_validated(b, check_finite=check_finite, as_inexact=True)
+
     # Validate shapes.
     if a1.shape[-1] != b1.shape[0]:
         raise ValueError("shapes of ab and b are not compatible.")
+
     (nlower, nupper) = l_and_u
     if nlower + nupper + 1 != a1.shape[0]:
         raise ValueError("invalid values for the number of lower and upper "
                          "diagonals: l+u+1 (%d) does not equal ab.shape[0] "
                          "(%d)" % (nlower + nupper + 1, ab.shape[0]))
 
+    # accommodate empty arrays
+    if b1.size == 0:
+        dt = solve(np.eye(1, dtype=a1.dtype), np.ones(1, dtype=b1.dtype)).dtype
+        return np.empty_like(b1, dtype=dt)
+
     overwrite_b = overwrite_b or _datacopied(b1, b)
     if a1.shape[-1] == 1:
         b2 = np.array(b1, copy=(not overwrite_b))
@@ -570,10 +588,16 @@ def solveh_banded(ab, b, overwrite_ab=False, overwrite_b=False, lower=False,
     """
     a1 = _asarray_validated(ab, check_finite=check_finite)
     b1 = _asarray_validated(b, check_finite=check_finite)
+
     # Validate shapes.
     if a1.shape[-1] != b1.shape[0]:
         raise ValueError("shapes of ab and b are not compatible.")
 
+    # accommodate empty arrays
+    if b1.size == 0:
+        dt = solve(np.eye(1, dtype=a1.dtype), np.ones(1, dtype=b1.dtype)).dtype
+        return np.empty_like(b1, dtype=dt)
+
     overwrite_b = overwrite_b or _datacopied(b1, b)
     overwrite_ab = overwrite_ab or _datacopied(a1, ab)
 
@@ -674,6 +698,10 @@ def solve_toeplitz(c_or_cr, b, check_finite=True):
     r, c, b, dtype, b_shape = _validate_args_for_toeplitz_ops(
         c_or_cr, b, check_finite, keep_b_shape=True)
 
+    # accommodate empty arrays
+    if b.size == 0:
+        return np.empty_like(b)
+
     # Form a 1-D array of values to be used in the matrix, containing a
     # reversed copy of r[1:], followed by c.
     vals = np.concatenate((r[-1:0:-1], c))
@@ -855,6 +883,10 @@ def solve_circulant(c, b, singular='raise', tol=None,
     if nc != nb:
         raise ValueError(f'Shapes of c {c.shape} and b {b.shape} are incompatible')
 
+    # accommodate empty arrays
+    if b.size == 0:
+        return np.empty_like(b)
+
     fc = np.fft.fft(np.moveaxis(c, caxis, -1), axis=-1)
     abs_fc = np.abs(fc)
     if tol is None:
@@ -940,6 +972,12 @@ def inv(a, overwrite_a=False, check_finite=True):
     a1 = _asarray_validated(a, check_finite=check_finite)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
         raise ValueError('expected square matrix')
+
+    # accommodate empty square matrices
+    if a1.size == 0:
+        dt = inv(np.eye(2, dtype=a1.dtype)).dtype
+        return np.empty_like(a1, dtype=dt)
+
     overwrite_a = overwrite_a or _datacopied(a1, a)
     getrf, getri, getri_lwork = get_lapack_funcs(('getrf', 'getri',
                                                   'getri_lwork'),
@@ -1049,9 +1087,9 @@ def det(a, overwrite_a=False, check_finite=True):
     # Empty array has determinant 1 because math.
     if min(*a1.shape) == 0:
         if a1.ndim == 2:
-            return np.float64(1.)
+            return a1.dtype.type(1.)
         else:
-            return np.ones(shape=a1.shape[:-2], dtype=np.float64)
+            return np.ones(shape=a1.shape[:-2], dtype=a1.dtype)
 
     # Scalar case
     if a1.shape[-2:] == (1, 1):
@@ -1422,7 +1460,7 @@ def pinv(a, *, atol=None, rtol=None, return_rank=False, check_finite=True,
     a = _asarray_validated(a, check_finite=check_finite)
     u, s, vh = _decomp_svd.svd(a, full_matrices=False, check_finite=False)
     t = u.dtype.char.lower()
-    maxS = np.max(s)
+    maxS = np.max(s, initial=0.)
 
     if rcond is not _NoValue or cond is not _NoValue:
         warn('Use of the "cond" and "rcond" keywords are deprecated and '
@@ -1525,7 +1563,7 @@ def pinvh(a, atol=None, rtol=None, lower=True, return_rank=False,
     a = _asarray_validated(a, check_finite=check_finite)
     s, u = _decomp.eigh(a, lower=lower, check_finite=False, driver='ev')
     t = u.dtype.char.lower()
-    maxS = np.max(np.abs(s))
+    maxS = np.max(np.abs(s), initial=0.)
 
     atol = 0. if atol is None else atol
     rtol = max(a.shape) * np.finfo(t).eps if (rtol is None) else rtol
@@ -1651,6 +1689,16 @@ def matrix_balance(A, permute=True, scale=True, separate=False,
     if not np.equal(*A.shape):
         raise ValueError('The data matrix for balancing should be square.')
 
+    # accommodate empty arrays
+    if A.size == 0:
+        b_n, t_n = matrix_balance(np.eye(2, dtype=A.dtype))
+        B = np.empty_like(A, dtype=b_n.dtype)
+        if separate:
+            scaling = np.ones_like(A, shape=len(A))
+            perm = np.arange(len(A))
+            return B, (scaling, perm)
+        return B, np.empty_like(A, dtype=t_n.dtype)
+
     gebal = get_lapack_funcs(('gebal'), (A,))
     B, lo, hi, ps, info = gebal(A, scale=scale, permute=permute,
                                 overwrite_a=overwrite_a)
@@ -1760,7 +1808,7 @@ def _validate_args_for_toeplitz_ops(c_or_cr, b, check_finite, keep_b_shape,
     if b.ndim == 1 and not keep_b_shape:
         b = b.reshape(-1, 1)
     elif b.ndim != 1:
-        b = b.reshape(b.shape[0], -1)
+        b = b.reshape(b.shape[0], -1 if b.size > 0 else 0)
 
     return r, c, b, dtype, b_shape
 
@@ -1898,6 +1946,11 @@ def matmul_toeplitz(c_or_cr, x, check_finite=False, workers=None):
     T_nrows = len(c)
     T_ncols = len(r)
     p = T_nrows + T_ncols - 1  # equivalent to len(embedded_col)
+    return_shape = (T_nrows,) if len(x_shape) == 1 else (T_nrows, m)
+
+    # accommodate empty arrays
+    if x.size == 0:
+        return np.empty_like(x, shape=return_shape)
 
     embedded_col = np.concatenate((c, r[-1:0:-1]))
 
@@ -1915,5 +1968,4 @@ def matmul_toeplitz(c_or_cr, x, check_finite=False, workers=None):
         mat_times_x = irfft(fft_mat*fft_x, axis=0,
                             workers=workers, n=p)[:T_nrows, :]
 
-    return_shape = (T_nrows,) if len(x_shape) == 1 else (T_nrows, m)
     return mat_times_x.reshape(*return_shape)
diff --git a/scipy/linalg/_decomp.py b/scipy/linalg/_decomp.py
index 659f0b88d4df..86f96db1c6a1 100644
--- a/scipy/linalg/_decomp.py
+++ b/scipy/linalg/_decomp.py
@@ -215,6 +215,22 @@ def eig(a, b=None, left=False, right=True, overwrite_a=False,
     a1 = _asarray_validated(a, check_finite=check_finite)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
         raise ValueError('expected square matrix')
+
+    # accommodate square empty matrices
+    if a1.size == 0:
+        w_n, vr_n = eig(numpy.eye(2, dtype=a1.dtype))
+        w = numpy.empty_like(a1, shape=(0,), dtype=w_n.dtype)
+        w = _make_eigvals(w, None, homogeneous_eigvals)
+        vl = numpy.empty_like(a1, shape=(0, 0), dtype=vr_n.dtype)
+        vr = numpy.empty_like(a1, shape=(0, 0), dtype=vr_n.dtype)
+        if not (left or right):
+            return w
+        if left:
+            if right:
+                return w, vl, vr
+            return w, vl
+        return w, vr
+
     overwrite_a = overwrite_a or (_datacopied(a1, a))
     if b is not None:
         b1 = _asarray_validated(b, check_finite=check_finite)
@@ -244,7 +260,6 @@ def eig(a, b=None, left=False, right=True, overwrite_a=False,
                                     compute_vl=compute_vl,
                                     compute_vr=compute_vr,
                                     overwrite_a=overwrite_a)
-        t = {'f': 'F', 'd': 'D'}[wr.dtype.char]
         w = wr + _I * wi
         w = _make_eigvals(w, None, homogeneous_eigvals)
 
@@ -469,6 +484,18 @@ def eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False,
     a1 = _asarray_validated(a, check_finite=check_finite)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
         raise ValueError('expected square "a" matrix')
+
+    # accommodate square empty matrices
+    if a1.size == 0:
+        w_n, v_n = eigh(numpy.eye(2, dtype=a1.dtype))
+
+        w = numpy.empty_like(a1, shape=(0,), dtype=w_n.dtype)
+        v = numpy.empty_like(a1, shape=(0, 0), dtype=v_n.dtype)
+        if eigvals_only:
+            return w
+        else:
+            return w, v
+
     overwrite_a = overwrite_a or (_datacopied(a1, a))
     cplx = True if iscomplexobj(a1) else False
     n = a1.shape[0]
@@ -592,8 +619,8 @@ def eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False,
             return w, v
     else:
         if info < -1:
-            raise LinAlgError('Illegal value in argument {} of internal {}'
-                              ''.format(-info, drv.typecode + pfx + driver))
+            raise LinAlgError(f'Illegal value in argument {-info} of internal '
+                              f'{drv.typecode + pfx + driver}')
         elif info > n:
             raise LinAlgError(f'The leading minor of order {info-n} of B is not '
                               'positive definite. The factorization of B '
@@ -781,8 +808,21 @@ def eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False,
 
     if len(a1.shape) != 2:
         raise ValueError('expected a 2-D array')
+
+    # accommodate square empty matrices
+    if a1.size == 0:
+        w_n, v_n = eig_banded(numpy.array([[0, 0], [1, 1]], dtype=a1.dtype))
+
+        w = numpy.empty_like(a1, shape=(0,), dtype=w_n.dtype)
+        v = numpy.empty_like(a1, shape=(0, 0), dtype=v_n.dtype)
+        if eigvals_only:
+            return w
+        else:
+            return w, v
+
     select, vl, vu, il, iu, max_ev = _check_select(
         select, select_range, max_ev, a1.shape[1])
+
     del select_range
     if select == 0:
         if a1.dtype.char in 'GFD':
@@ -1447,6 +1487,17 @@ def hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True):
         raise ValueError('expected square matrix')
     overwrite_a = overwrite_a or (_datacopied(a1, a))
 
+    if a1.size == 0:
+        h3 = hessenberg(numpy.eye(3, dtype=a1.dtype))
+        h = numpy.empty(a1.shape, dtype=h3.dtype)
+        if not calc_q:
+            return h
+        else:
+            h3, q3 = hessenberg(numpy.eye(3, dtype=a1.dtype), calc_q=True)
+            q = numpy.empty(a1.shape, dtype=q3.dtype)
+            h = numpy.empty(a1.shape, dtype=h3.dtype)
+            return h, q
+
     # if 2x2 or smaller: already in Hessenberg
     if a1.shape[0] <= 2:
         if calc_q:
diff --git a/scipy/linalg/_decomp_cholesky.py b/scipy/linalg/_decomp_cholesky.py
index 76de9de3fab7..9bd1101f18e0 100644
--- a/scipy/linalg/_decomp_cholesky.py
+++ b/scipy/linalg/_decomp_cholesky.py
@@ -1,6 +1,7 @@
 """Cholesky decomposition functions."""
 
-from numpy import asarray_chkfinite, asarray, atleast_2d
+import numpy as np
+from numpy import asarray_chkfinite, asarray, atleast_2d, empty_like
 
 # Local imports
 from ._misc import LinAlgError, _datacopied
@@ -27,7 +28,8 @@ def _cholesky(a, lower=False, overwrite_a=False, clean=True,
 
     # Quick return for square empty array
     if a1.size == 0:
-        return a1.copy(), lower
+        dt = cholesky(np.eye(1, dtype=a1.dtype)).dtype
+        return empty_like(a1, dtype=dt), lower
 
     overwrite_a = overwrite_a or _datacopied(a1, a)
     potrf, = get_lapack_funcs(('potrf',), (a1,))
@@ -198,11 +200,18 @@ def cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True):
     else:
         b1 = asarray(b)
         c = asarray(c)
+
     if c.ndim != 2 or c.shape[0] != c.shape[1]:
         raise ValueError("The factored matrix c is not square.")
     if c.shape[1] != b1.shape[0]:
         raise ValueError(f"incompatible dimensions ({c.shape} and {b1.shape})")
 
+    # accommodate empty arrays
+    if b1.size == 0:
+        dt = cho_solve((np.eye(2, dtype=b1.dtype), True),
+                        np.ones(2, dtype=c.dtype)).dtype
+        return empty_like(b1, dtype=dt)
+
     overwrite_b = overwrite_b or _datacopied(b1, b)
 
     potrs, = get_lapack_funcs(('potrs',), (c, b1))
@@ -278,6 +287,11 @@ def cholesky_banded(ab, overwrite_ab=False, lower=False, check_finite=True):
     else:
         ab = asarray(ab)
 
+    # accommodate square empty matrices
+    if ab.size == 0:
+        dt = cholesky_banded(np.array([[0, 0], [1, 1]], dtype=ab.dtype)).dtype
+        return empty_like(ab, dtype=dt)
+
     pbtrf, = get_lapack_funcs(('pbtrf',), (ab,))
     c, info = pbtrf(ab, lower=lower, overwrite_ab=overwrite_ab)
     if info > 0:
@@ -346,6 +360,12 @@ def cho_solve_banded(cb_and_lower, b, overwrite_b=False, check_finite=True):
     if cb.shape[-1] != b.shape[0]:
         raise ValueError("shapes of cb and b are not compatible.")
 
+    # accommodate empty arrays
+    if b.size == 0:
+        m = cholesky_banded(np.array([[0, 0], [1, 1]], dtype=cb.dtype))
+        dt = cho_solve_banded((m, True), np.ones(2, dtype=b.dtype)).dtype
+        return empty_like(b, dtype=dt)
+
     pbtrs, = get_lapack_funcs(('pbtrs',), (cb, b))
     x, info = pbtrs(cb, b, lower=lower, overwrite_b=overwrite_b)
     if info > 0:
diff --git a/scipy/linalg/_decomp_lu.py b/scipy/linalg/_decomp_lu.py
index 2bb088b45317..f6926eead031 100644
--- a/scipy/linalg/_decomp_lu.py
+++ b/scipy/linalg/_decomp_lu.py
@@ -106,7 +106,15 @@ def lu_factor(a, overwrite_a=False, check_finite=True):
         a1 = asarray_chkfinite(a)
     else:
         a1 = asarray(a)
+
+    # accommodate empty arrays
+    if a1.size == 0:
+        lu = np.empty_like(a1)
+        piv = np.arange(0, dtype=np.int32)
+        return lu, piv
+
     overwrite_a = overwrite_a or (_datacopied(a1, a))
+
     getrf, = get_lapack_funcs(('getrf',), (a1,))
     lu, piv, info = getrf(a1, overwrite_a=overwrite_a)
     if info < 0:
@@ -171,10 +179,17 @@ def lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True):
         b1 = asarray_chkfinite(b)
     else:
         b1 = asarray(b)
+
     overwrite_b = overwrite_b or _datacopied(b1, b)
+
     if lu.shape[0] != b1.shape[0]:
         raise ValueError(f"Shapes of lu {lu.shape} and b {b1.shape} are incompatible")
 
+    # accommodate empty arrays
+    if b1.size == 0:
+        m = lu_solve((np.eye(2, dtype=lu.dtype), [0, 1]), np.ones(2, dtype=b.dtype))
+        return np.empty_like(b1, dtype=m.dtype)
+
     getrs, = get_lapack_funcs(('getrs',), (lu, b1))
     x, info = getrs(lu, piv, b1, trans=trans, overwrite_b=overwrite_b)
     if info == 0:
diff --git a/scipy/linalg/_decomp_qr.py b/scipy/linalg/_decomp_qr.py
index 2b158015e7bf..942ce1b3f8b5 100644
--- a/scipy/linalg/_decomp_qr.py
+++ b/scipy/linalg/_decomp_qr.py
@@ -131,7 +131,34 @@ def qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False,
         a1 = numpy.asarray(a)
     if len(a1.shape) != 2:
         raise ValueError("expected a 2-D array")
+
     M, N = a1.shape
+
+    # accommodate empty arrays
+    if a1.size == 0:
+        K = min(M, N)
+
+        if mode not in ['economic', 'raw']:
+            Q = numpy.empty_like(a1, shape=(M, M))
+            Q[...] = numpy.identity(M)
+            R = numpy.empty_like(a1)
+        else:
+            Q = numpy.empty_like(a1, shape=(M, K))
+            R = numpy.empty_like(a1, shape=(K, N))
+
+        if pivoting:
+            Rj = R, numpy.arange(N, dtype=numpy.int32)
+        else:
+            Rj = R,
+
+        if mode == 'r':
+            return Rj
+        elif mode == 'raw':
+            qr = numpy.empty_like(a1, shape=(M, N))
+            tau = numpy.zeros_like(a1, shape=(K,))
+            return ((qr, tau),) + Rj
+        return (Q,) + Rj
+
     overwrite_a = overwrite_a or (_datacopied(a1, a))
 
     if pivoting:
@@ -279,6 +306,10 @@ def qr_multiply(a, c, mode='right', pivoting=False, conjugate=False,
     raw = qr(a, overwrite_a, None, "raw", pivoting)
     Q, tau = raw[0]
 
+    # accommodate empty arrays
+    if c.size == 0:
+        return (numpy.empty_like(c),) + raw[1:]
+
     gor_un_mqr, = get_lapack_funcs(('ormqr',), (Q,))
     if gor_un_mqr.typecode in ('s', 'd'):
         trans = "T"
@@ -399,7 +430,25 @@ def rq(a, overwrite_a=False, lwork=None, mode='full', check_finite=True):
         a1 = numpy.asarray(a)
     if len(a1.shape) != 2:
         raise ValueError('expected matrix')
+
     M, N = a1.shape
+
+    # accommodate empty arrays
+    if a1.size == 0:
+        K = min(M, N)
+
+        if not mode == 'economic':
+            R = numpy.empty_like(a1)
+            Q = numpy.empty_like(a1, shape=(N, N))
+            Q[...] = numpy.identity(N)
+        else:
+            R = numpy.empty_like(a1, shape=(M, K))
+            Q = numpy.empty_like(a1, shape=(K, N))
+
+        if mode == 'r':
+            return R
+        return R, Q
+
     overwrite_a = overwrite_a or (_datacopied(a1, a))
 
     gerqf, = get_lapack_funcs(('gerqf',), (a1,))
diff --git a/scipy/linalg/_decomp_schur.py b/scipy/linalg/_decomp_schur.py
index ed41e4f9ea7e..32a9190040d2 100644
--- a/scipy/linalg/_decomp_schur.py
+++ b/scipy/linalg/_decomp_schur.py
@@ -124,14 +124,24 @@ def schur(a, output='real', lwork=None, overwrite_a=False, sort=None,
         a1 = asarray(a, dtype=numpy.dtype("long"))
     if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
         raise ValueError('expected square matrix')
+
     typ = a1.dtype.char
     if output in ['complex', 'c'] and typ not in ['F', 'D']:
         if typ in _double_precision:
             a1 = a1.astype('D')
-            typ = 'D'
         else:
             a1 = a1.astype('F')
-            typ = 'F'
+
+    # accommodate empty matrix
+    if a1.size == 0:
+        t0, z0 = schur(numpy.eye(2, dtype=a1.dtype))
+        if sort is None:
+            return (numpy.empty_like(a1, dtype=t0.dtype),
+                    numpy.empty_like(a1, dtype=z0.dtype))
+        else:
+            return (numpy.empty_like(a1, dtype=t0.dtype),
+                    numpy.empty_like(a1, dtype=z0.dtype), 0)
+
     overwrite_a = overwrite_a or (_datacopied(a1, a))
     gees, = get_lapack_funcs(('gees',), (a1,))
     if lwork is None or lwork == -1:
@@ -176,7 +186,7 @@ def sfunction(x):
     elif info > 0:
         raise LinAlgError("Schur form not found. Possibly ill-conditioned.")
 
-    if sort_t == 0:
+    if sort is None:
         return result[0], result[-3]
     else:
         return result[0], result[-3], result[1]
diff --git a/scipy/linalg/_decomp_svd.py b/scipy/linalg/_decomp_svd.py
index a2452e3dc21b..1902f60710af 100644
--- a/scipy/linalg/_decomp_svd.py
+++ b/scipy/linalg/_decomp_svd.py
@@ -109,6 +109,25 @@ def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False,
     if len(a1.shape) != 2:
         raise ValueError('expected matrix')
     m, n = a1.shape
+
+    # accommodate empty matrix
+    if a1.size == 0:
+        u0, s0, v0 = svd(numpy.eye(2, dtype=a1.dtype))
+
+        s = numpy.empty_like(a1, shape=(0,), dtype=s0.dtype)
+        if full_matrices:
+            u = numpy.empty_like(a1, shape=(m, m), dtype=u0.dtype)
+            u[...] = numpy.identity(m)
+            v = numpy.empty_like(a1, shape=(n, n), dtype=v0.dtype)
+            v[...] = numpy.identity(n)
+        else:
+            u = numpy.empty_like(a1, shape=(m, 0), dtype=u0.dtype)
+            v = numpy.empty_like(a1, shape=(0, n), dtype=v0.dtype)
+        if compute_uv:
+            return u, s, v
+        else:
+            return s
+
     overwrite_a = overwrite_a or (_datacopied(a1, a))
 
     if not isinstance(lapack_driver, str):
@@ -236,14 +255,8 @@ def svdvals(a, overwrite_a=False, check_finite=True):
     array([ 1.,  1.,  1.,  1.])
 
     """
-    a = _asarray_validated(a, check_finite=check_finite)
-    if a.size:
-        return svd(a, compute_uv=0, overwrite_a=overwrite_a,
-                   check_finite=False)
-    elif len(a.shape) != 2:
-        raise ValueError('expected matrix')
-    else:
-        return numpy.empty(0)
+    return svd(a, compute_uv=0, overwrite_a=overwrite_a,
+               check_finite=check_finite)
 
 
 def diagsvd(s, M, N):
@@ -340,7 +353,7 @@ def orth(A, rcond=None):
     M, N = u.shape[0], vh.shape[1]
     if rcond is None:
         rcond = numpy.finfo(s.dtype).eps * max(M, N)
-    tol = numpy.amax(s) * rcond
+    tol = numpy.amax(s, initial=0.) * rcond
     num = numpy.sum(s > tol, dtype=int)
     Q = u[:, :num]
     return Q
@@ -404,7 +417,7 @@ def null_space(A, rcond=None):
     M, N = u.shape[0], vh.shape[1]
     if rcond is None:
         rcond = numpy.finfo(s.dtype).eps * max(M, N)
-    tol = numpy.amax(s) * rcond
+    tol = numpy.amax(s, initial=0.) * rcond
     num = numpy.sum(s > tol, dtype=int)
     Q = vh[num:,:].T.conj()
     return Q
diff --git a/scipy/linalg/_matfuncs.py b/scipy/linalg/_matfuncs.py
index 64fc65df291d..aaad628c90ad 100644
--- a/scipy/linalg/_matfuncs.py
+++ b/scipy/linalg/_matfuncs.py
@@ -856,6 +856,12 @@ def khatri_rao(a, b):
         raise ValueError("The number of columns for both arrays "
                          "should be equal.")
 
+    # accommodate empty arrays
+    if a.size == 0 or b.size == 0:
+        m = a.shape[0] * b.shape[0]
+        n = a.shape[1]
+        return np.empty_like(a, shape=(m, n))
+
     # c = np.vstack([np.kron(a[:, k], b[:, k]) for k in range(b.shape[1])]).T
     c = a[..., :, np.newaxis, :] * b[..., np.newaxis, :, :]
     return c.reshape((-1,) + c.shape[2:])
diff --git a/scipy/linalg/_matfuncs_inv_ssq.py b/scipy/linalg/_matfuncs_inv_ssq.py
index d84f53dc6c83..1decffae2e52 100644
--- a/scipy/linalg/_matfuncs_inv_ssq.py
+++ b/scipy/linalg/_matfuncs_inv_ssq.py
@@ -368,7 +368,7 @@ def _inverse_squaring_helper(T0, theta):
     tmp_diag = np.diag(T)
     if np.count_nonzero(tmp_diag) != n:
         raise Exception('Diagonal entries of T must be nonzero')
-    while np.max(np.absolute(tmp_diag - 1)) > theta[7]:
+    while np.max(np.absolute(tmp_diag - 1), initial=0.) > theta[7]:
         tmp_diag = np.sqrt(tmp_diag)
         s0 += 1
 
@@ -761,7 +761,7 @@ def _logm_triu(T):
     # Construct T0 with the appropriate type,
     # depending on the dtype and the spectrum of T.
     T_diag = np.diag(T)
-    keep_it_real = np.isrealobj(T) and np.min(T_diag) >= 0
+    keep_it_real = np.isrealobj(T) and np.min(T_diag, initial=0.) >= 0
     if keep_it_real:
         T0 = T
     else:
@@ -866,7 +866,7 @@ def _logm(A):
     try:
         if np.array_equal(A, np.triu(A)):
             A = _logm_force_nonsingular_triangular_matrix(A)
-            if np.min(np.diag(A)) < 0:
+            if np.min(np.diag(A), initial=0.) < 0:
                 A = A.astype(complex)
             return _logm_triu(A)
         else:
diff --git a/scipy/linalg/_matfuncs_sqrtm.py b/scipy/linalg/_matfuncs_sqrtm.py
index 9c84ca573bd2..e76c08932eb4 100644
--- a/scipy/linalg/_matfuncs_sqrtm.py
+++ b/scipy/linalg/_matfuncs_sqrtm.py
@@ -51,7 +51,7 @@ def _sqrtm_triu(T, blocksize=64):
 
     """
     T_diag = np.diag(T)
-    keep_it_real = np.isrealobj(T) and np.min(T_diag) >= 0
+    keep_it_real = np.isrealobj(T) and np.min(T_diag, initial=0.) >= 0
 
     # Cast to complex as necessary + ensure double precision
     if not keep_it_real:
diff --git a/scipy/linalg/_special_matrices.py b/scipy/linalg/_special_matrices.py
index df9b767481d1..791ef4a6ba5d 100644
--- a/scipy/linalg/_special_matrices.py
+++ b/scipy/linalg/_special_matrices.py
@@ -336,6 +336,12 @@ def kron(a, b):
            [3, 3, 3, 4, 4, 4]])
 
     """
+    # accommodate empty arrays
+    if a.size == 0 or b.size == 0:
+        m = a.shape[0] * b.shape[0]
+        n = a.shape[1] * b.shape[1]
+        return np.empty_like(a, shape=(m, n))
+
     if not a.flags['CONTIGUOUS']:
         a = np.reshape(a, a.shape)
     if not b.flags['CONTIGUOUS']:
diff --git a/scipy/linalg/tests/meson.build b/scipy/linalg/tests/meson.build
index cc0728ba84a4..deae9788a977 100644
--- a/scipy/linalg/tests/meson.build
+++ b/scipy/linalg/tests/meson.build
@@ -17,7 +17,6 @@ python_sources = [
   'test_lapack.py',
   'test_matfuncs.py',
   'test_matmul_toeplitz.py',
-  'test_misc.py',
   'test_procrustes.py',
   'test_sketches.py',
   'test_solve_toeplitz.py',
diff --git a/scipy/linalg/tests/test_basic.py b/scipy/linalg/tests/test_basic.py
index 67610035d059..f127a6ff1e50 100644
--- a/scipy/linalg/tests/test_basic.py
+++ b/scipy/linalg/tests/test_basic.py
@@ -175,6 +175,23 @@ def test_native_list_arguments(self):
         x = solve_banded((l, u), ab, b)
         assert_array_almost_equal(dot(a, x), b)
 
+    @pytest.mark.parametrize('dt_ab', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt_ab, dt_b):
+        # ab contains one empty row corresponding to the diagonal
+        ab = np.array([[]], dtype=dt_ab)
+        b = np.array([], dtype=dt_b)
+        x = solve_banded((0, 0), ab, b)
+
+        assert x.shape == (0,)
+        assert x.dtype == solve(np.eye(1, dtype=dt_ab), np.ones(1, dtype=dt_b)).dtype
+
+        b = np.empty((0, 0), dtype=dt_b)
+        x = solve_banded((0, 0), ab, b)
+
+        assert x.shape == (0, 0)
+        assert x.dtype == solve(np.eye(1, dtype=dt_ab), np.ones(1, dtype=dt_b)).dtype
+
 
 class TestSolveHBanded:
 
@@ -490,6 +507,22 @@ def test_native_list_arguments(self):
         x = solveh_banded(ab, b)
         assert_array_almost_equal(x, [0.0, 1.0, 0.0, 0.0])
 
+    @pytest.mark.parametrize('dt_ab', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt_ab, dt_b):
+        # ab contains one empty row corresponding to the diagonal
+        ab = np.array([[]], dtype=dt_ab)
+        b = np.array([], dtype=dt_b)
+        x = solveh_banded(ab, b)
+
+        assert x.shape == (0,)
+        assert x.dtype == solve(np.eye(1, dtype=dt_ab), np.ones(1, dtype=dt_b)).dtype
+
+        b = np.empty((0, 0), dtype=dt_b)
+        x = solveh_banded(ab, b)
+
+        assert x.shape == (0, 0)
+        assert x.dtype == solve(np.eye(1, dtype=dt_ab), np.ones(1, dtype=dt_b)).dtype
 
 class TestSolve:
     def setup_method(self):
@@ -737,13 +770,6 @@ def test_ill_condition_warning(self):
             warnings.simplefilter('error')
             assert_raises(LinAlgWarning, solve, a, b)
 
-    def test_empty_rhs(self):
-        a = np.eye(2)
-        b = [[], []]
-        x = solve(a, b)
-        assert_(x.size == 0, 'Returned array is not empty')
-        assert_(x.shape == (2, 0), 'Returned empty array shape is wrong')
-
     def test_multiple_rhs(self):
         a = np.eye(2)
         b = np.random.rand(2, 3, 4)
@@ -821,6 +847,24 @@ def test_all_type_size_routine_combinations(self):
                                 rtol=tol * size,
                                 err_msg=err_msg)
 
+    @pytest.mark.parametrize('dt_a', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt_a, dt_b):
+        a = np.empty((0, 0), dtype=dt_a)
+        b = np.empty(0, dtype=dt_b)
+        x = solve(a, b)
+
+        assert x.size == 0
+        dt_nonempty = solve(np.eye(2, dtype=dt_a), np.ones(2, dtype=dt_b)).dtype
+        assert x.dtype == dt_nonempty
+
+    def test_empty_rhs(self):
+        a = np.eye(2)
+        b = [[], []]
+        x = solve(a, b)
+        assert_(x.size == 0, 'Returned array is not empty')
+        assert_(x.shape == (2, 0), 'Returned empty array shape is wrong')
+
 
 class TestSolveTriangular:
 
@@ -883,6 +927,26 @@ def test_check_finite(self):
         sol = solve_triangular(A, b, lower=True, check_finite=False)
         assert_array_almost_equal(sol, [1, 0])
 
+    @pytest.mark.parametrize('dt_a', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt_a, dt_b):
+        a = np.empty((0, 0), dtype=dt_a)
+        b = np.empty(0, dtype=dt_b)
+        x = solve_triangular(a, b)
+
+        assert x.size == 0
+        dt_nonempty = solve_triangular(
+            np.eye(2, dtype=dt_a), np.ones(2, dtype=dt_b)
+        ).dtype
+        assert x.dtype == dt_nonempty
+
+    def test_empty_rhs(self):
+        a = np.eye(2)
+        b = [[], []]
+        x = solve_triangular(a, b)
+        assert_(x.size == 0, 'Returned array is not empty')
+        assert_(x.shape == (2, 0), 'Returned empty array shape is wrong')
+
 
 class TestInv:
     def setup_method(self):
@@ -926,6 +990,13 @@ def test_check_finite(self):
         a_inv = inv(a, check_finite=False)
         assert_array_almost_equal(dot(a, a_inv), [[1, 0], [0, 1]])
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        a_inv = inv(a)
+        assert a_inv.size == 0
+        assert a_inv.dtype == inv(np.eye(2, dtype=dt)).dtype
+
 
 class TestDet:
     def setup_method(self):
@@ -1043,6 +1114,14 @@ def test_empty_edge_cases(self):
         with assert_raises(ValueError, match='Last 2 dimensions'):
             det(np.array([[[]]]))
 
+    @pytest.mark.xfail(reason="det(float32) -> float64, is it a bug?")
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty_dtype(self, dt):
+        a = np.empty((3, 0, 0), dtype=dt)
+        d = det(a)
+        assert d.shape == (3,)
+        assert d.dtype == det(np.empty((3, 1, 1), dtype=dt)).dtype
+
     def test_overwrite_a(self):
         # If all conditions are met then input should be overwritten;
         #   - dtype is one of 'fdFD'
@@ -1348,7 +1427,7 @@ def test_check_finite(self):
                             atol=25 * _eps_cast(a.dtype),
                             err_msg="driver: %s" % lapack_driver)
 
-    def test_zero_size(self):
+    def test_empty(self):
         for a_shape, b_shape in (((0, 2), (0,)),
                                  ((0, 4), (0, 2)),
                                  ((4, 0), (4,)),
@@ -1362,6 +1441,16 @@ def test_zero_size(self):
             assert_(rank == 0, 'expected rank 0')
             assert_equal(s, np.empty((0,)))
 
+    @pytest.mark.parametrize('dt_a', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty_dtype(self, dt_a, dt_b):
+        a = np.empty((0, 0), dtype=dt_a)
+        b = np.empty(0, dtype=dt_b)
+        x, residues, rank, s = lstsq(a, b)
+
+        assert x.size == 0
+        dt_nonempty = lstsq(np.eye(2, dtype=dt_a), np.ones(2, dtype=dt_b))[0].dtype
+        assert x.dtype == dt_nonempty
 
 class TestPinv:
     def setup_method(self):
@@ -1457,6 +1546,13 @@ def test_positional_deprecation(self):
         with pytest.deprecated_call(match="use keyword arguments"):
             pinv(np.ones((2,2)), 0., 1e-10)
 
+    @pytest.mark.parametrize('dt', [float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        a_pinv = pinv(a)
+        assert a_pinv.size == 0
+        assert a_pinv.dtype == pinv(np.eye(2, dtype=dt)).dtype
+
 
 class TestPinvSymmetric:
 
@@ -1528,6 +1624,13 @@ def test_atol_rtol(self):
         assert_allclose(norm(adiff1), 1e-4, rtol=0.1)
         assert_allclose(norm(adiff2), 1e-4, rtol=0.1)
 
+    @pytest.mark.parametrize('dt', [float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        a_pinv = pinvh(a)
+        assert a_pinv.size == 0
+        assert a_pinv.dtype == pinv(np.eye(2, dtype=dt)).dtype
+
 
 @pytest.mark.parametrize('scale', (1e-20, 1., 1e20))
 @pytest.mark.parametrize('pinv_', (pinv, pinvh))
@@ -1644,6 +1747,17 @@ def test_keepdims_kwd(self):
         assert_allclose(b, c)
         assert_(b.shape == c.shape)
 
+    def test_empty(self):
+        a = np.empty((0, 0))
+        assert_allclose(norm(a), 0.)
+        assert_allclose(norm(a, axis=0), np.zeros((0,)))
+        assert_allclose(norm(a, keepdims=True), np.zeros((1, 1)))
+
+        a = np.empty((0, 3))
+        assert_allclose(norm(a), 0.)
+        assert_allclose(norm(a, axis=0), np.zeros((3,)))
+        assert_allclose(norm(a, keepdims=True), np.zeros((1, 1)))
+
 
 class TestOverwrite:
     def test_solve(self):
@@ -1759,6 +1873,16 @@ def test_native_list_arguments(self):
         y = solve(circulant(c), b)
         assert_allclose(x, y)
 
+    def test_empty(self):
+        c = []
+        b = []
+        x = solve_circulant(c, b)
+        assert_allclose(x, [])
+
+        b = np.empty((0, 0))
+        x = solve_circulant(c, b)
+        assert_allclose(x, np.empty((0, 0)))
+
 
 class TestMatrix_Balance:
 
@@ -1823,3 +1947,27 @@ def test_perm_and_scaling(self):
             ip[p] = np.arange(A.shape[0])
             assert_allclose(y, np.diag(s)[ip, :])
             assert_allclose(solve(y, A).dot(y), x)
+
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        b, t = matrix_balance(a)
+
+        assert b.size == 0
+        assert t.size == 0
+
+        b_n, t_n = matrix_balance(np.eye(2, dtype=dt))
+        assert b.dtype == b_n.dtype
+        assert t.dtype == t_n.dtype
+
+
+        b, (scale, perm) = matrix_balance(a, separate=True)
+        assert b.size == 0
+        assert scale.size == 0
+        assert perm.size == 0
+
+        b_n, (scale_n, perm_n) = matrix_balance(a, separate=True)
+        assert b.dtype == b_n.dtype
+        assert scale.dtype == scale_n.dtype
+        assert perm.dtype == perm_n.dtype
+
diff --git a/scipy/linalg/tests/test_decomp.py b/scipy/linalg/tests/test_decomp.py
index f9e347050627..b5a350d82eb4 100644
--- a/scipy/linalg/tests/test_decomp.py
+++ b/scipy/linalg/tests/test_decomp.py
@@ -116,6 +116,17 @@ def test_finite(self):
         exact_w = [(9+sqrt(93))/2, 0, (9-sqrt(93))/2]
         assert_array_almost_equal(w, exact_w)
 
+    @pytest.mark.parametrize('dt', [int, float, float32, complex, complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        w = eigvals(a)
+        assert w.shape == (0,)
+        assert w.dtype == eigvals(np.eye(2, dtype=dt)).dtype
+
+        w = eigvals(a, homogeneous_eigvals=True)
+        assert w.shape == (2, 0)
+        assert w.dtype == eigvals(np.eye(2, dtype=dt)).dtype
+
 
 class TestEig:
 
@@ -368,6 +379,28 @@ def test_gh_11577(self):
             assert np.isclose(D, 4.0, atol=1e-14).any()
             assert np.isclose(D, 8.0, atol=1e-14).any()
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        w, vr = eig(a)
+
+        w_n, vr_n = eig(np.eye(2, dtype=dt))
+
+        assert w.shape == (0,)
+        assert w.dtype == w_n.dtype  #eigvals(np.eye(2, dtype=dt)).dtype
+
+        assert_allclose(vr, np.empty((0, 0)))
+        assert vr.shape == (0, 0)
+        assert vr.dtype == vr_n.dtype
+
+        w, vr = eig(a, homogeneous_eigvals=True)
+        assert w.shape == (2, 0)
+        assert w.dtype == w_n.dtype
+
+        assert vr.shape == (0, 0)
+        assert vr.dtype == vr_n.dtype
+
+
 
 class TestEigBanded:
     def setup_method(self):
@@ -648,6 +681,22 @@ def test_zgbtrs(self):
         y_lin = linalg.solve(self.comp_mat, self.bc)
         assert_array_almost_equal(y, y_lin)
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a_band = np.empty((0, 0), dtype=dt)
+        w, v = eig_banded(a_band)
+
+        w_n, v_n = eig_banded(np.array([[0, 0], [1, 1]], dtype=dt))
+
+        assert w.shape == (0,)
+        assert w.dtype == w_n.dtype
+
+        assert v.shape == (0, 0)
+        assert v.dtype == v_n.dtype
+
+        w = eig_banded(a_band, eigvals_only=True)
+        assert w.shape == (0,)
+        assert w.dtype == w_n.dtype
 
 class TestEigTridiagonal:
     def setup_method(self):
@@ -936,6 +985,25 @@ def test_deprecation_results(self):
         assert_allclose(v_dep, v)
 
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        w, v = eigh(a)
+
+        w_n, v_n = eigh(np.eye(2, dtype=dt))
+
+        assert w.shape == (0,)
+        assert w.dtype == w_n.dtype
+
+        assert v.shape == (0, 0)
+        assert v.dtype == v_n.dtype
+
+        w = eigh(a, eigvals_only=True)
+        assert_allclose(w, np.empty((0,)))
+
+        assert w.shape == (0,)
+        assert w.dtype == w_n.dtype
+
 class TestSVD_GESDD:
     lapack_driver = 'gesdd'
 
@@ -1084,6 +1152,64 @@ def test_large_matrix(self):
         assert_allclose(s[0], 1.0)
         assert_allclose(u[0, 0] * vh[0, -1], 1.0)
 
+    @pytest.mark.parametrize("m", [0, 1, 2])
+    @pytest.mark.parametrize("n", [0, 1, 2])
+    @pytest.mark.parametrize('dtype', DTYPES)
+    def test_shape_dtype(self, m, n, dtype):
+        a = np.zeros((m, n), dtype=dtype)
+        k = min(m, n)
+        dchar = a.dtype.char
+        real_dchar = dchar.lower() if dchar in 'FD' else dchar
+
+        u, s, v = svd(a)
+        assert_equal(u.shape, (m, m))
+        assert_equal(u.dtype, dtype)
+        assert_equal(s.shape, (k,))
+        assert_equal(s.dtype, np.dtype(real_dchar))
+        assert_equal(v.shape, (n, n))
+        assert_equal(v.dtype, dtype)
+
+        u, s, v = svd(a, full_matrices=False)
+        assert_equal(u.shape, (m, k))
+        assert_equal(u.dtype, dtype)
+        assert_equal(s.shape, (k,))
+        assert_equal(s.dtype, np.dtype(real_dchar))
+        assert_equal(v.shape, (k, n))
+        assert_equal(v.dtype, dtype)
+
+        s = svd(a, compute_uv=False)
+        assert_equal(s.shape, (k,))
+        assert_equal(s.dtype, np.dtype(real_dchar))
+
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize(("m", "n"), [(0, 0), (0, 2), (2, 0)])
+    def test_empty(self, dt, m, n):
+        a0 = np.eye(3, dtype=dt)
+        u0, s0, v0 = svd(a0)
+
+        a = np.empty((m, n), dtype=dt)
+        u, s, v = svd(a)
+        assert_allclose(u, np.identity(m))
+        assert_allclose(s, np.empty((0,)))
+        assert_allclose(v, np.identity(n))
+
+        assert u.dtype == u0.dtype
+        assert v.dtype == v0.dtype
+        assert s.dtype == s0.dtype
+
+        u, s, v = svd(a, full_matrices=False)
+        assert_allclose(u, np.empty((m, 0)))
+        assert_allclose(s, np.empty((0,)))
+        assert_allclose(v, np.empty((0, n)))
+
+        assert u.dtype == u0.dtype
+        assert v.dtype == v0.dtype
+        assert s.dtype == s0.dtype
+
+        s = svd(a, compute_uv=False)
+        assert_allclose(s, np.empty((0,)))
+
+        assert s.dtype == s0.dtype
 
 class TestSVD_GESVD(TestSVD_GESDD):
     lapack_driver = 'gesvd'
@@ -1099,11 +1225,16 @@ def test_svd_gesdd_nofegfault():
 
 class TestSVDVals:
 
-    def test_empty(self):
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
         for a in [[]], np.empty((2, 0)), np.ones((0, 3)):
+            a = np.array(a, dtype=dt)
             s = svdvals(a)
             assert_equal(s, np.empty(0))
 
+            s0 = svdvals(np.eye(2, dtype=dt))
+            assert s.dtype == s0.dtype
+
     def test_simple(self):
         a = [[1, 2, 3], [1, 2, 3], [2, 5, 6]]
         s = svdvals(a)
@@ -1714,6 +1845,100 @@ def test_lwork(self):
         assert_raises(Exception, qr, (a,), {'lwork': 0})
         assert_raises(Exception, qr, (a,), {'lwork': 2})
 
+    @pytest.mark.parametrize("m", [0, 1, 2])
+    @pytest.mark.parametrize("n", [0, 1, 2])
+    @pytest.mark.parametrize("pivoting", [False, True])
+    @pytest.mark.parametrize('dtype', DTYPES)
+    def test_shape_dtype(self, m, n, pivoting, dtype):
+        k = min(m, n)
+
+        a = np.zeros((m, n), dtype=dtype)
+        q, r, *other = qr(a, pivoting=pivoting)
+        assert_equal(q.shape, (m, m))
+        assert_equal(q.dtype, dtype)
+        assert_equal(r.shape, (m, n))
+        assert_equal(r.dtype, dtype)
+        assert len(other) == (1 if pivoting else 0)
+        if pivoting:
+            p, = other
+            assert_equal(p.shape, (n,))
+            assert_equal(p.dtype, np.int32)
+
+        r, *other = qr(a, mode='r', pivoting=pivoting)
+        assert_equal(r.shape, (m, n))
+        assert_equal(r.dtype, dtype)
+        assert len(other) == (1 if pivoting else 0)
+        if pivoting:
+            p, = other
+            assert_equal(p.shape, (n,))
+            assert_equal(p.dtype, np.int32)
+
+        q, r, *other = qr(a, mode='economic', pivoting=pivoting)
+        assert_equal(q.shape, (m, k))
+        assert_equal(q.dtype, dtype)
+        assert_equal(r.shape, (k, n))
+        assert_equal(r.dtype, dtype)
+        assert len(other) == (1 if pivoting else 0)
+        if pivoting:
+            p, = other
+            assert_equal(p.shape, (n,))
+            assert_equal(p.dtype, np.int32)
+
+        (raw, tau), r, *other = qr(a, mode='raw', pivoting=pivoting)
+        assert_equal(raw.shape, (m, n))
+        assert_equal(raw.dtype, dtype)
+        assert_equal(tau.shape, (k,))
+        assert_equal(tau.dtype, dtype)
+        assert_equal(r.shape, (k, n))
+        assert_equal(r.dtype, dtype)
+        assert len(other) == (1 if pivoting else 0)
+        if pivoting:
+            p, = other
+            assert_equal(p.shape, (n,))
+            assert_equal(p.dtype, np.int32)
+
+    @pytest.mark.parametrize(("m", "n"), [(0, 0), (0, 2), (2, 0)])
+    def test_empty(self, m, n):
+        k = min(m, n)
+
+        a = np.empty((m, n))
+        q, r = qr(a)
+        assert_allclose(q, np.identity(m))
+        assert_allclose(r, np.empty((m, n)))
+
+        q, r, p = qr(a, pivoting=True)
+        assert_allclose(q, np.identity(m))
+        assert_allclose(r, np.empty((m, n)))
+        assert_allclose(p, np.arange(n))
+
+        r, = qr(a, mode='r')
+        assert_allclose(r, np.empty((m, n)))
+
+        q, r = qr(a, mode='economic')
+        assert_allclose(q, np.empty((m, k)))
+        assert_allclose(r, np.empty((k, n)))
+
+        (raw, tau), r = qr(a, mode='raw')
+        assert_allclose(raw, np.empty((m, n)))
+        assert_allclose(tau, np.empty((k,)))
+        assert_allclose(r, np.empty((k, n)))
+
+    def test_multiply_empty(self):
+        a = np.empty((0, 0))
+        c = np.empty((0, 0))
+        cq, r = qr_multiply(a, c)
+        assert_allclose(cq, np.empty((0, 0)))
+
+        a = np.empty((0, 2))
+        c = np.empty((2, 0))
+        cq, r = qr_multiply(a, c)
+        assert_allclose(cq, np.empty((2, 0)))
+
+        a = np.empty((2, 0))
+        c = np.empty((0, 2))
+        cq, r = qr_multiply(a, c)
+        assert_allclose(cq, np.empty((0, 2)))
+
 
 class TestRQ:
     def test_simple(self):
@@ -1820,6 +2045,45 @@ def test_check_finite(self):
         assert_array_almost_equal(q @ q.T, eye(3))
         assert_array_almost_equal(r @ q, a)
 
+    @pytest.mark.parametrize("m", [0, 1, 2])
+    @pytest.mark.parametrize("n", [0, 1, 2])
+    @pytest.mark.parametrize('dtype', DTYPES)
+    def test_shape_dtype(self, m, n, dtype):
+        k = min(m, n)
+
+        a = np.zeros((m, n), dtype=dtype)
+        r, q = rq(a)
+        assert_equal(q.shape, (n, n))
+        assert_equal(r.shape, (m, n))
+        assert_equal(r.dtype, dtype)
+        assert_equal(q.dtype, dtype)
+
+        r = rq(a, mode='r')
+        assert_equal(r.shape, (m, n))
+        assert_equal(r.dtype, dtype)
+
+        r, q = rq(a, mode='economic')
+        assert_equal(r.shape, (m, k))
+        assert_equal(r.dtype, dtype)
+        assert_equal(q.shape, (k, n))
+        assert_equal(q.dtype, dtype)
+
+    @pytest.mark.parametrize(("m", "n"), [(0, 0), (0, 2), (2, 0)])
+    def test_empty(self, m, n):
+        k = min(m, n)
+
+        a = np.empty((m, n))
+        r, q = rq(a)
+        assert_allclose(r, np.empty((m, n)))
+        assert_allclose(q, np.identity(n))
+
+        r = rq(a, mode='r')
+        assert_allclose(r, np.empty((m, n)))
+
+        r, q = rq(a, mode='economic')
+        assert_allclose(r, np.empty((m, k)))
+        assert_allclose(q, np.empty((k, n)))
+
 
 class TestSchur:
 
@@ -1875,6 +2139,23 @@ def test_check_finite(self):
         t, z = schur(a, check_finite=False)
         assert_array_almost_equal(z @ t @ z.conj().T, a)
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        t, z = schur(a)
+        t0, z0 = schur(np.eye(2, dtype=dt))
+        assert_allclose(t, np.empty((0, 0)))
+        assert_allclose(z, np.empty((0, 0)))
+        assert t.dtype == t0.dtype
+        assert z.dtype == z0.dtype
+
+        t, z, sdim = schur(a, sort='lhp')
+        assert_allclose(t, np.empty((0, 0)))
+        assert_allclose(z, np.empty((0, 0)))
+        assert_equal(sdim, 0)
+        assert t.dtype == t0.dtype
+        assert z.dtype == z0.dtype
+
 
 class TestHessenberg:
 
@@ -1952,6 +2233,21 @@ def test_2x2(self):
         assert_array_almost_equal(q2, np.eye(2))
         assert_array_almost_equal(h2, b)
 
+    @pytest.mark.parametrize('dt', [int, float, float32, complex, complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        h = hessenberg(a)
+        assert h.shape == (0, 0)
+        assert h.dtype == hessenberg(np.eye(3, dtype=dt)).dtype
+
+        h, q = hessenberg(a, calc_q=True)
+        h3, q3 = hessenberg(a, calc_q=True)
+        assert h.shape == (0, 0)
+        assert h.dtype == h3.dtype
+
+        assert q.shape == (0, 0)
+        assert q.dtype == q3.dtype
+
 
 blas_provider = blas_version = None
 if CONFIG is not None:
@@ -2596,6 +2892,15 @@ def test_orth():
     for dt, n in itertools.product(dtypes, sizes):
         _check_orth(n, dt)
 
+@pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+def test_orth_empty(dt):
+    a = np.empty((0, 0), dtype=dt)
+    a0 = np.eye(2, dtype=dt)
+
+    oa = orth(a)
+    assert oa.dtype == orth(a0).dtype
+    assert oa.shape == (0, 0)
+
 
 def test_null_space():
     np.random.seed(1)
@@ -2635,6 +2940,16 @@ def test_null_space():
             assert_equal(Y.shape, (n, n - 6))
 
 
+@pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+def test_null_space_empty(dt):
+    a = np.empty((0, 0), dtype=dt)
+    a0 = np.eye(2, dtype=dt)
+    nsa = null_space(a)
+
+    assert nsa.shape == (0, 0)
+    assert nsa.dtype == null_space(a0).dtype
+
+
 def test_subspace_angles():
     H = hadamard(8, float)
     A = H[:, :3]
@@ -2692,6 +3007,17 @@ def test_subspace_angles():
     assert_allclose(subspace_angles(a, b), 0., atol=1e-14)
     assert_allclose(subspace_angles(b, a), 0., atol=1e-14)
 
+    # Empty
+    a = np.empty((0, 0))
+    b = np.empty((0, 0))
+    assert_allclose(subspace_angles(a, b), np.empty((0,)))
+    a = np.empty((2, 0))
+    b = np.empty((2, 0))
+    assert_allclose(subspace_angles(a, b), np.empty((0,)))
+    a = np.empty((0, 2))
+    b = np.empty((0, 3))
+    assert_allclose(subspace_angles(a, b), np.empty((0,)))
+
 
 class TestCDF2RDF:
 
diff --git a/scipy/linalg/tests/test_decomp_cholesky.py b/scipy/linalg/tests/test_decomp_cholesky.py
index 9354bf93a967..61bbc7e544f1 100644
--- a/scipy/linalg/tests/test_decomp_cholesky.py
+++ b/scipy/linalg/tests/test_decomp_cholesky.py
@@ -1,12 +1,12 @@
 import pytest
-from numpy.testing import assert_array_almost_equal, assert_array_equal
+import numpy as np
+from numpy.testing import assert_array_almost_equal
 from pytest import raises as assert_raises
 
-import numpy as np
 from numpy import array, transpose, dot, conjugate, zeros_like, empty
 from numpy.random import random
-from scipy.linalg import cholesky, cholesky_banded, cho_solve_banded, \
-     cho_factor, cho_solve
+from scipy.linalg import (cholesky, cholesky_banded, cho_solve_banded,
+     cho_factor, cho_solve)
 
 from scipy.linalg._testutils import assert_no_overwrite
 
@@ -82,6 +82,35 @@ def test_int_overflow(self):
 
        cholesky(x, check_finite=False, overwrite_a=True)  # should not segfault
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt, dt_b):
+        a = empty((0, 0), dtype=dt)
+
+        c = cholesky(a)
+        assert c.shape == (0, 0)
+        assert c.dtype == cholesky(np.eye(2, dtype=dt)).dtype
+
+        c_and_lower = (c, True)
+        b = np.asarray([], dtype=dt_b)
+        x = cho_solve(c_and_lower, b)
+        assert x.shape == (0,)
+        assert x.dtype == cho_solve((np.eye(2, dtype=dt), True),
+                                     np.ones(2, dtype=dt_b)).dtype
+
+        b = empty((0, 0), dtype=dt_b)
+        x = cho_solve(c_and_lower, b)
+        assert x.shape == (0, 0)
+        assert x.dtype == cho_solve((np.eye(2, dtype=dt), True),
+                                     np.ones(2, dtype=dt_b)).dtype
+
+        a1 = array([])
+        a2 = array([[]])
+        a3 = []
+        a4 = [[]]
+        for x in ([a1, a2, a3, a4]):
+            assert_raises(ValueError, cholesky, x)
+
 
 class TestCholeskyBanded:
     """Tests for cholesky_banded() and cho_solve_banded."""
@@ -181,6 +210,30 @@ def test_lower_complex(self):
         x = cho_solve_banded((c, True), b)
         assert_array_almost_equal(x, [0.0, 0.0, 1.0j, 1.0])
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt, dt_b):
+        ab = empty((0, 0), dtype=dt)
+
+        cb = cholesky_banded(ab)
+        assert cb.shape == (0, 0)
+
+        m = cholesky_banded(np.array([[0, 0], [1, 1]], dtype=dt))
+        assert cb.dtype == m.dtype
+
+        cb_and_lower = (cb, True)
+        b = np.asarray([], dtype=dt_b)
+        x = cho_solve_banded(cb_and_lower, b)
+        assert x.shape == (0,)
+
+        dtype_nonempty = cho_solve_banded((m, True), np.ones(2, dtype=dt_b)).dtype
+        assert x.dtype == dtype_nonempty
+
+        b = empty((0, 0), dtype=dt_b)
+        x = cho_solve_banded(cb_and_lower, b)
+        assert x.shape == (0, 0)
+        assert x.dtype == dtype_nonempty
+
 
 class TestOverwrite:
     def test_cholesky(self):
@@ -203,17 +256,13 @@ def test_cho_solve_banded(self):
         assert_no_overwrite(lambda b: cho_solve_banded((xcho, False), b),
                             [(3,)])
 
+class TestChoFactor:
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        x, lower = cho_factor(a)
 
-class TestEmptyArray:
-    def test_cho_factor_empty_square(self):
-        a = empty((0, 0))
-        b = array([])
-        c = array([[]])
-        d = []
-        e = [[]]
-
-        x, _ = cho_factor(a)
-        assert_array_equal(x, a)
+        assert x.shape == (0, 0)
 
-        for x in ([b, c, d, e]):
-            assert_raises(ValueError, cho_factor, x)
+        xx, lower = cho_factor(np.eye(2, dtype=dt))
+        assert x.dtype == xx.dtype
diff --git a/scipy/linalg/tests/test_decomp_lu.py b/scipy/linalg/tests/test_decomp_lu.py
index d17ef62f5e3b..b94ef7ffb1c5 100644
--- a/scipy/linalg/tests/test_decomp_lu.py
+++ b/scipy/linalg/tests/test_decomp_lu.py
@@ -3,7 +3,12 @@
 
 import numpy as np
 from scipy.linalg import lu, lu_factor, lu_solve, get_lapack_funcs, solve
-from numpy.testing import assert_allclose, assert_array_equal
+from numpy.testing import assert_allclose, assert_array_equal, assert_equal
+
+
+REAL_DTYPES = [np.float32, np.float64]
+COMPLEX_DTYPES = [np.complex64, np.complex128]
+DTYPES = REAL_DTYPES + COMPLEX_DTYPES
 
 
 class TestLU:
@@ -242,6 +247,26 @@ def test_simple_known(self):
             assert_allclose(LU, np.array([[2, 1], [0, 1]]))
             assert_array_equal(P, np.array([0, 1]))
 
+    @pytest.mark.parametrize("m", [0, 1, 2])
+    @pytest.mark.parametrize("n", [0, 1, 2])
+    @pytest.mark.parametrize('dtype', DTYPES)
+    def test_shape_dtype(self, m, n,  dtype):
+        k = min(m, n)
+
+        a = np.eye(m, n, dtype=dtype)
+        lu, p = lu_factor(a)
+        assert_equal(lu.shape, (m, n))
+        assert_equal(lu.dtype, dtype)
+        assert_equal(p.shape, (k,))
+        assert_equal(p.dtype, np.int32)
+
+    @pytest.mark.parametrize(("m", "n"), [(0, 0), (0, 2), (2, 0)])
+    def test_empty(self, m, n):
+        a = np.zeros((m, n))
+        lu, p = lu_factor(a)
+        assert_allclose(lu, np.empty((m, n)))
+        assert_allclose(p, np.arange(0))
+
 
 class TestLUSolve:
     def setup_method(self):
@@ -265,3 +290,19 @@ def test_check_finite(self):
         lu_a = lu_factor(a, check_finite=False)
         x2 = lu_solve(lu_a, b, check_finite=False)
         assert_allclose(x1, x2)
+
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    @pytest.mark.parametrize('dt_b', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt, dt_b):
+        lu_and_piv = (np.empty((0, 0), dtype=dt), np.array([]))
+        b = np.asarray([], dtype=dt_b)
+        x = lu_solve(lu_and_piv, b)
+        assert x.shape == (0,)
+
+        m = lu_solve((np.eye(2, dtype=dt), [0, 1]), np.ones(2, dtype=dt_b))
+        assert x.dtype == m.dtype
+
+        b = np.empty((0, 0), dtype=dt_b)
+        x = lu_solve(lu_and_piv, b)
+        assert x.shape == (0, 0)
+        assert x.dtype == m.dtype
diff --git a/scipy/linalg/tests/test_decomp_polar.py b/scipy/linalg/tests/test_decomp_polar.py
index 44c4c9c09ed6..b5dc75255288 100644
--- a/scipy/linalg/tests/test_decomp_polar.py
+++ b/scipy/linalg/tests/test_decomp_polar.py
@@ -29,6 +29,9 @@
     [[1, 2, 3+2j], [3, 4-1j, -4j]],
     [[1, 2], [3-2j, 4+0.5j], [5, 5]],
     [[10000, 10, 1], [-1, 2, 3j], [0, 1, 2]],
+    np.empty((0, 0)),
+    np.empty((0, 2)),
+    np.empty((2, 0)),
 ]
 
 
diff --git a/scipy/linalg/tests/test_matfuncs.py b/scipy/linalg/tests/test_matfuncs.py
index b50122f81fa8..89e290478415 100644
--- a/scipy/linalg/tests/test_matfuncs.py
+++ b/scipy/linalg/tests/test_matfuncs.py
@@ -238,6 +238,17 @@ def test_readonly(self):
         a.flags.writeable = False
         logm(a)
 
+    @pytest.mark.xfail(reason="ValueError: attempt to get argmax of an empty sequence")
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        log_a = logm(a)
+        a0 = np.eye(2, dtype=dt)
+        log_a0 = logm(a0)
+
+        assert log_a.shape == (0, 0)
+        assert log_a.dtype == log_a0.dtype
+
 
 class TestSqrtM:
     def test_round_trip_random_float(self):
@@ -490,6 +501,16 @@ def test_data_size_preservation_comp_in_comp_out(self):
             M = np.array([[2j, 4], [0, -2j]], dtype=np.complex256)
             assert sqrtm(M).dtype == np.complex256
 
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty(self, dt):
+        a = np.empty((0, 0), dtype=dt)
+        s = sqrtm(a)
+        a0 = np.eye(2, dtype=dt)
+        s0 = sqrtm(a0)
+
+        assert s.shape == (0, 0)
+        assert s.dtype == s0.dtype
+
 
 class TestFractionalMatrixPower:
     def test_round_trip_random_complex(self):
@@ -685,11 +706,13 @@ def test_single_elt(self):
         elt = expm(1)
         assert_allclose(elt, np.array([[np.e]]))
 
-    def test_empty_matrix_input(self):
+    @pytest.mark.parametrize('dt', [int, float, np.float32, complex, np.complex64])
+    def test_empty_matrix_input(self, dt):
         # handle gh-11082
-        A = np.zeros((0, 0))
+        A = np.zeros((0, 0), dtype=dt)
         result = expm(A)
         assert result.size == 0
+        assert result.dtype == A.dtype
 
     def test_2x2_input(self):
         E = np.e
@@ -1011,3 +1034,14 @@ def test_equality_of_two_equations(self):
                           for k in range(b.shape[1])]).T
 
         assert_array_equal(res1, res2)
+
+    def test_empty(self):
+        a = np.empty((0, 2))
+        b = np.empty((3, 2))
+        res = khatri_rao(a, b)
+        assert_allclose(res, np.empty((0, 2)))
+
+        a = np.empty((3, 0))
+        b = np.empty((5, 0))
+        res = khatri_rao(a, b)
+        assert_allclose(res, np.empty((15, 0)))
diff --git a/scipy/linalg/tests/test_matmul_toeplitz.py b/scipy/linalg/tests/test_matmul_toeplitz.py
index b480e9d398ba..609c0043c099 100644
--- a/scipy/linalg/tests/test_matmul_toeplitz.py
+++ b/scipy/linalg/tests/test_matmul_toeplitz.py
@@ -94,6 +94,15 @@ def test_complex(self):
             self.rng.normal(size=(n//2, 3))*1j
         self.do(x, c, r, False)
 
+    def test_empty(self):
+        c = []
+        r = []
+        x = []
+        self.do(x, c, r, False)
+
+        x = np.empty((0, 0))
+        self.do(x, c, r, False)
+
     def test_exceptions(self):
 
         n = 100
diff --git a/scipy/linalg/tests/test_misc.py b/scipy/linalg/tests/test_misc.py
deleted file mode 100644
index 1c10923e088b..000000000000
--- a/scipy/linalg/tests/test_misc.py
+++ /dev/null
@@ -1,5 +0,0 @@
-from scipy.linalg import norm
-
-
-def test_norm():
-    assert norm([]) == 0.0
diff --git a/scipy/linalg/tests/test_procrustes.py b/scipy/linalg/tests/test_procrustes.py
index f41fd0e254c7..bb5049bbdb42 100644
--- a/scipy/linalg/tests/test_procrustes.py
+++ b/scipy/linalg/tests/test_procrustes.py
@@ -189,3 +189,13 @@ def test_orthogonal_procrustes_skbio_example():
     B_approx = scale * np.dot(A, R) + B_mu
     assert_allclose(B_approx, B_orig)
     assert_allclose(B / norm(B), B_standardized)
+
+
+def test_empty():
+    a = np.empty((0, 0))
+    r, s = orthogonal_procrustes(a, a)
+    assert_allclose(r, np.empty((0, 0)))
+
+    a = np.empty((0, 3))
+    r, s = orthogonal_procrustes(a, a)
+    assert_allclose(r, np.identity(3))
diff --git a/scipy/linalg/tests/test_solve_toeplitz.py b/scipy/linalg/tests/test_solve_toeplitz.py
index ecced19e2d39..ebae76f35c9d 100644
--- a/scipy/linalg/tests/test_solve_toeplitz.py
+++ b/scipy/linalg/tests/test_solve_toeplitz.py
@@ -119,3 +119,14 @@ def test_unstable():
 
     assert_allclose(solution1, solution2)
 
+
+def test_empty():
+    c = []
+    r = []
+    b = []
+    x = solve_toeplitz((c, r), b)
+    assert_allclose(x, [])
+
+    b = np.empty((0, 0))
+    x = solve_toeplitz((c, r), b)
+    assert_allclose(x, np.empty((0, 0)))
diff --git a/scipy/linalg/tests/test_special_matrices.py b/scipy/linalg/tests/test_special_matrices.py
index 3edc8c176717..dd36ddf0ccba 100644
--- a/scipy/linalg/tests/test_special_matrices.py
+++ b/scipy/linalg/tests/test_special_matrices.py
@@ -219,6 +219,12 @@ def test_basic(self):
                           [33, 44]])
         assert_array_equal(a, expected)
 
+    def test_empty(self):
+        m1 = np.empty((0, 2))
+        m2 = np.empty((1, 3))
+        a = kron(m1, m2)
+        assert_allclose(a, np.empty((0, 6)))
+
 
 class TestHelmert:
 


Code review comments:
diff_hunk: @@ -1043,6 +1114,14 @@ def test_empty_edge_cases(self):
         with assert_raises(ValueError, match='Last 2 dimensions'):
             det(np.array([[[]]]))
 
+    @pytest.mark.xfail(reason="det(float32) -> float64, is it a bug?")
path: scipy/linalg/tests/test_basic.py
text: @user1:
This looks like it's on purpose:
```python
    # ... avoid overflows for single precision
    # Cf. np.linalg.det(np.diag([1e+38, 1e+38]).astype(np.float32))
    dtype_char = a1.dtype.char
    if dtype_char in 'fF':
        dtype_char = 'd' if dtype_char.islower() else 'D'
```

The documented return type is `float or complex` - not super clear, but it's on purpose. `det` upcasting to 64-bit isn't unreasonable I think. So I suggest improving the docs a bit and updating this test so it passes.

@author:
Good catch, thanks! 
Indeed, it's even documented, https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.det.html

>    Even the input array is single precision (float32 or complex64), the result
    will be returned in double precision (float64 or complex128) to prevent
    overflows.

PR updated to make it consistent for empty and non-empty arrays.

diff_hunk: @@ -1049,9 +1087,9 @@ def det(a, overwrite_a=False, check_finite=True):
     # Empty array has determinant 1 because math.
     if min(*a1.shape) == 0:
         if a1.ndim == 2:
-            return np.float64(1.)
+            return a1.dtype.type(1.)
path: scipy/linalg/_basic.py
text: @user1:
Related to the `xfail` comment: this change and the one two lines below should probably be reverted, since `float64` was intended.

diff_hunk: @@ -855,6 +883,10 @@ def solve_circulant(c, b, singular='raise', tol=None,
     if nc != nb:
         raise ValueError(f'Shapes of c {c.shape} and b {b.shape} are incompatible')
 
+    # accommodate empty arrays
+    if b.size == 0:
+        return np.empty_like(b)
path: scipy/linalg/_basic.py
text: @user1:
minor issue: this misses the dtype determination. I commented on the corresponding test case with more detail.

diff_hunk: @@ -1759,6 +1873,16 @@ def test_native_list_arguments(self):
         y = solve(circulant(c), b)
         assert_allclose(x, y)
 
+    def test_empty(self):
+        c = []
+        b = []
+        x = solve_circulant(c, b)
+        assert_allclose(x, [])
+
+        b = np.empty((0, 0))
+        x = solve_circulant(c, b)
+        assert_allclose(x, np.empty((0, 0)))
path: scipy/linalg/tests/test_basic.py
text: @user1:
A small issue here:
```python
>>> c = np.array([2, 2, 4])
>>> linalg.solve_circulant(c, np.empty((3,0), dtype=np.int64))
array([], shape=(3, 0), dtype=int64)
>>> linalg.solve_circulant(c, np.empty((3,1), dtype=np.int64)).dtype
dtype('float64')
```

Output should always be real or complex floating-point

@author:
Updated, thanks!