FIX: Fix known bugs in MarkovChain (markov/mc_tools.jl) - #392
Conversation
- gth_solve: throw DimensionMismatch on non-square input, instead of silently returning an incorrect result - gth_solve!: avoid temporary arrays from row/column slicing and out-of-place normalization (at n=200: 1201 allocations down to 5) - check_stochastic_matrix: replace the fixed 5e-15 absolute row-sum tolerance, which rejects valid large matrices, with a size-dependent one; this also makes Float32 matrices usable - test_mc_tools.jl: remove a type-pirated Base.isapprox overload that never dispatched due to type invariance (and would recurse infinitely if it did) Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
There was a problem hiding this comment.
Pull request overview
This PR addresses correctness and performance issues in QuantEcon.jl’s Markov chain utilities (src/markov/mc_tools.jl), focused on making gth_solve/gth_solve! safer for invalid inputs and making stochastic-matrix validation scale appropriately with matrix size and element type.
Changes:
- Add an explicit squareness check to
gth_solve!(and thereforegth_solve) to throwDimensionMismatchinstead of producing incorrect/NaN-contaminated results. - Reduce allocations in
gth_solve!by switching to views and in-place scaling/normalization. - Replace the fixed row-sum tolerance in
check_stochastic_matrixwith a size- and eltype-dependent tolerance, and update tests accordingly (also removing a non-functional/type-piratingBase.isapproxoverload in tests).
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated no comments.
| File | Description |
|---|---|
src/markov/mc_tools.jl |
Adds dimension validation to gth_solve!, removes avoidable allocations via views/in-place ops, and scales stochastic-matrix row-sum tolerance with n * eps(T) (with a small-matrix floor). |
test/test_mc_tools.jl |
Removes an invalid Base.isapprox overload and adds tests for non-square gth_solve inputs plus the updated row-sum tolerance behavior (including a Float32 case). |
Corrections to the check_stochastic_matrix tolerance from review: - Cap the tolerance at sqrt(eps(T)) so that it can never approach the scale of the entries: n * eps(Float16) reaches exactly 1 at n=1024, which accepted an all-zero Float16 matrix. - For sparse matrices, scale with the maximum number of stored entries per row instead of n: the rounding error of a row does not grow with the total number of states, and n * eps(Float32) at n=100_000 accepted rows 1% short of 1. - Accumulate half/single precision row sums in Float64, so that large legitimately normalized low-precision matrices are not rejected by rounding introduced by the check itself. - Fall back to eps(Float64) for abstract element types, for which eps(T) raises a MethodError. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
|
A review by ChatGPT caught three genuine flaws in the tolerance introduced here, fixed in 16f3575. (1) 🤖 Generated with Claude Code (Claude Fable 5) |
This PR fixes several bugs in
markov/mc_tools.jlfound in a code review, following up on #385/#391.Fixes
gth_solve/gth_solve!silently returned an incorrect result (NaN-contaminated) for non-square input, since onlysize(A, 1)was read and the elimination runs under@inbounds. ADimensionMismatchis now thrown.gth_solve!allocated temporaries in the elimination loop (row-slice copy in the pivot-scale sum and an out-of-place column scaling) and in the final normalization. These are now views/in-place: at n=200, allocations drop from 1201 (842 KiB) to 5 (322 KiB, essentially the input copy ingth_solve), with a ~10% speedup at small sizes (benchmarkSUITE["mc_tools"]["gth_solve"], Add benchmarks for MarkovChain (markov/mc_tools.jl) #391).check_stochastic_matrixused a fixed absolute tolerance of 5e-15 on row sums, which rejects valid large matrices: both the rounding of the entries (e.g. from row normalization) and the accumulation in the recomputed sum grow linearly with n, so e.g. a row-normalized random 1000x1000 matrix typically fails the check. The tolerance is nowmax(5e-15, n * eps(T))forT<:AbstractFloat(witheps(Float64)otherwise); the floor keeps the current behavior for small matrices. This also makesFloat32transition matrices usable, which the fixed tolerance rejected essentially always.test_mc_tools.jldefined aBase.isapproxoverload forVector{Vector{<:Real}}— type piracy that in fact never dispatched (by type invariance,Vector{Vector{Float64}} <: Vector{Vector{<:Real}}does not hold), and whose body compared the outer arguments, so it would have recursed infinitely if it ever ran. It is removed; the tests exercise the genericisapproxfallback as they always effectively did.Tests
New test sets cover the non-square
DimensionMismatch(for theFloat64,Int, and in-place entry points), row-sum errors within and beyond the size-dependent tolerance, and construction from a row-normalizedFloat32matrix.🤖 Generated with Claude Code (Claude Fable 5)