Jonathan Wenger

We need `from_svd` as well for `LinearOperators` where one can compute the SVD but not the cholesky decomposition.

Further, always keep the sampling dimensions as the first dimensions, such that a statement like `np.shape(sample_)` has the value of `size` as its first components.

- [ ] Write tests to test all of the above (already for distributions or just random variables?)

> I would second this. One question: Does this only concern `SymmetricMatrixVariateNormal`, or do we also support the non-symmetric version? (just curious, not saying we need to support the non-symmetric...

@fxbriol Could you link the website with the test functions for quadrature, that you mentioned offline?

I agree with your assessment @marvinpfoertner. This behaviour is particularly nasty ```python >>> x = np.random.randn(2,1) >>> n.pdf(x) array([0.04105293, 0.11377797]) ``` I guess there is some implicit NumPy broadcasting that's...

The intended behavior for any methods (`.pdf`, `.cdf`, etc.) is that they take arguments of shape `(n_0, n_1, ...., d_0, d_1, ...)` where `n_i` determines the batches and `d_i` the...

# To Do - [ ] Fix broadcasting bug as described above - [ ] Adjust docstring to include the information about input shape - [ ] Write doctests for...

Having `BlockDiagonal` as a `LinearOperator` seems very convenient but is there a difference in efficiency or is this just an alias for the implemented `Kronecker(A=Identity(num_copies), B=block)`?

It's very surprising that it is actually slower than the dense form of NumPy. That looks like a bug in the `Kronecker` implementation to me.