You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
It seems that Diagonal matrices in LinearAlgebra do not fully implement the interface of pinv. There is an implementation of pinv(D::Diagonal)here and of pinv(D::Diagonal, rtol) just below. There is no implementation of a pseudo-inverse with regularization, as in pinv(D::Diagonal; atol = ..., rtol = ...).
On the other hand, the generic routine for dense arrays explicitly checks for diagonal matrices and implements a custom pinv for that case with thresholds here. However, this code does not return a diagonal matrix:
Might anyone else be interested in having a more complete implementation of pinv for Diagonal matrices which returns a Diagonal?
Combining the implementation in dense.jl with the one in diagonal.jl, it could look something like this:
function LinearAlgebra.pinv(D::Diagonal{T}; atol::Real=0.0, rtol::Real= (eps(real(float(oneunit(T))))*length(D.diag))*iszero(atol)) where T
Di =similar(D.diag, typeof(inv(oneunit(T))))
if!isempty(D.diag)
maxabsA =maximum(abs, D.diag)
abstol =max(rtol * maxabsA, atol)
for i in1:length(D.diag)
ifabs(D.diag[i]) > abstol
invD =inv(D.diag[i])
ifisfinite(invD)
Di[i] = invD
continueendend# fallback
Di[i] =zero(T)
endendDiagonal(Di)
end
There has been quite a discussion about default choices of the regularization parameters (e.g. #652, #444 and issues referenced there). The implementation above does not change that. It is merely about returning a Diagonal matrix to get this behaviour:
On the same topic: for consistency it may be worth adding atol and rtol keyword arguments to pinv(x::Number) as well, returning the same as it would for the corresponding 1-by-1 matrix.
It seems that Diagonal matrices in LinearAlgebra do not fully implement the interface of
pinv
. There is an implementation ofpinv(D::Diagonal)
here and ofpinv(D::Diagonal, rtol)
just below. There is no implementation of a pseudo-inverse with regularization, as inpinv(D::Diagonal; atol = ..., rtol = ...)
.On the other hand, the generic routine for dense arrays explicitly checks for diagonal matrices and implements a custom pinv for that case with thresholds here. However, this code does not return a diagonal matrix:
Might anyone else be interested in having a more complete implementation of pinv for Diagonal matrices which returns a Diagonal?
Combining the implementation in dense.jl with the one in diagonal.jl, it could look something like this:
There has been quite a discussion about default choices of the regularization parameters (e.g. #652, #444 and issues referenced there). The implementation above does not change that. It is merely about returning a Diagonal matrix to get this behaviour:
The text was updated successfully, but these errors were encountered: