Use exact algorithm for cond(A,1), and fix for singular matrices (J…

…uliaLang#33548)

NEWS.md

@@ -58,6 +58,8 @@ Standard library changes
 
 * `ldlt` and non-pivoted `lu` now throw a new `ZeroPivotException` type ([#33372]).
 
+* `cond(A, p)` with `p=1` or `p=Inf` now computes the exact condition number instead of an estimate ([#33547]).
+
 #### Random
 
 * `AbstractRNG`s now behave like scalars when used in broadcasting ([#33213]).

stdlib/LinearAlgebra/src/dense.jl

@@ -1418,15 +1418,22 @@ function cond(A::AbstractMatrix, p::Real=2)
  if p == 2
  v = svdvals(A)
  maxv = maximum(v)
- return maxv == 0.0 ? oftype(real(A[1,1]),Inf) : maxv / minimum(v)
+ return iszero(maxv) ? oftype(real(maxv), Inf) : maxv / minimum(v)
  elseif p == 1 || p == Inf
  checksquare(A)
- return _cond1Inf(A, p)
+ try
+ Ainv = inv(A)
+ return opnorm(A, p)*opnorm(Ainv, p)
+ catch e
+ if isa(e, LAPACKException) || isa(e, SingularException)
+ return convert(float(real(eltype(A))), Inf)
+ else
+ rethrow()
+ end
+ end
  end
  throw(ArgumentError("p-norm must be 1, 2 or Inf, got $p"))
 end
-_cond1Inf(A::StridedMatrix{<:BlasFloat}, p::Real) = _cond1Inf(lu(A), p, opnorm(A, p))
-_cond1Inf(A::AbstractMatrix, p::Real) = opnorm(A, p)*opnorm(inv(A), p)
 
 ## Lyapunov and Sylvester equation
 

stdlib/LinearAlgebra/src/lu.jl

@@ -478,13 +478,6 @@ inv!(A::LU{T,<:StridedMatrix}) where {T} =
  ldiv!(A.factors, copy(A), Matrix{T}(I, size(A, 1), size(A, 1)))
 inv(A::LU{<:BlasFloat,<:StridedMatrix}) = inv!(copy(A))
 
-function _cond1Inf(A::LU{<:BlasFloat,<:StridedMatrix}, p::Number, normA::Real)
- if p != 1 && p != Inf
- throw(ArgumentError("p must be either 1 or Inf"))
- end
- return inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, normA))
-end
-
 # Tridiagonal
 
 # See dgttrf.f

stdlib/LinearAlgebra/test/dense.jl

@@ -29,6 +29,16 @@ Random.seed!(1234321)
  @test_throws ArgumentError cond(a,3)
  end
  end
+ @testset "Singular matrices" for p in (1, 2, Inf)
+ @test cond(zeros(Int, 2, 2), p) == Inf
+ @test cond(zeros(2, 2), p) == Inf
+ @test cond([0 0; 1 1], p) == Inf
+ @test cond([0. 0.; 1. 1.], p) == Inf
+ end
+ @testset "Issue #33547, condition number of 2x2 matrix" begin
+ M = [1.0 -2.0; -2.0 -1.5]
+ @test cond(M, 1) ≈ 2.227272727272727
+ end
 end
 
 areal = randn(n,n)/2

stdlib/LinearAlgebra/test/lu.jl

@@ -265,10 +265,6 @@ end
  @test_throws DomainError logdet([1 1; 1 -1])
 end
 
-@testset "Issue 21453" begin
- @test_throws ArgumentError LinearAlgebra._cond1Inf(lu(randn(5,5)), 2, 2.0)
-end
-
 @testset "REPL printing" begin
  bf = IOBuffer()
  show(bf, "text/plain", lu(Matrix(I, 4, 4)))