Skip to content

Levenberg-Marquardt with Sparse J: J isn't square?! #599

New issue
New issue
Closed
Closed
Levenberg-Marquardt with Sparse J: J isn't square?!#599
JuliaArrays/ArrayInterface.jl
#462
 (+1)
Labels
bugSomething isn't working
@luke-kiernan

Description

@luke-kiernan
luke-kiernan
opened on May 6, 2025 · edited by luke-kiernan

Describe the bug 🐞

I'm trying to run the LevenbergMarquardt solver on the system of equations arising from a power flow problem. After setting everything up, NonlinearSolve.jl is saying that the Jacobian matrix isn't square. All matrices I've provided are indeed square.

Minimal Reproducible Example 👇

using NonlinearSolve
using SciMLBase
using SparseArrays
using PowerFlows, PowerSystemCaseBuilder
const PF = PowerFlows
const PSB = PowerSystemCaseBuilder
# small system for demo purposes: actual target system has 10k buses. 
sys = PSB.build_system(PSB.MatpowerTestSystems, "matpower_case5_sys")
pf = PF.ACPowerFlow()
data = PF.PowerFlowData(pf, sys)
residual = PF.ACPowerFlowResidual(data, 1)
J = PF.ACPowerFlowJacobian(data, 1) # this initializes the sparse structure.
J0 = deepcopy(J.Jv) # take a copy of the sparse structure.
x0 = PF.calculate_x0(data, 1)
@assert size(J0, 1) == size(J0, 2) # yes, dimensions match. 
@assert size(J0, 1) == size(x0, 1)

function f!(du::Vector{Float64}, u::Vector{Float64}, ::Any)
    println(size(du))
    println(size(u))
    residual(du, u, 1) # in-place compute f at u
end

function jac!(A::SparseMatrixCSC{Float64, Int32}, u::Vector{Float64}, ::Any)
    println(size(u))
    println(size(A))
    PF.update_data!(data, u, 1)
    J(A, 1) # in-place compute J at the point stored in data (set to u above)
end

func = NonlinearFunction{true, SciMLBase.AutoSpecialize}(
    f!;
    jac = jac!,
    jac_prototype = J0,
)

prob = NonlinearProblem(func, x0)
solver = NonlinearSolveFirstOrder.LevenbergMarquardt()
solve(prob, solver)

Error & Stacktrace ⚠️

ERROR: LoadError: ArgumentError: KLU only accepts square matrices.
Stacktrace:
  [1] LinearSolveSparseArraysExt.KLU.KLUFactorization(A::SparseMatrixCSC{Float64, Int32})
    @ LinearSolveSparseArraysExt.KLU ~/.julia/packages/LinearSolve/O6GSm/src/KLU/klu.jl:177
  [2] init_cacheval
    @ ~/.julia/packages/LinearSolve/O6GSm/ext/LinearSolveSparseArraysExt.jl:145 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:308 [inlined]
  [4] init_cacheval(alg::LinearSolve.DefaultLinearSolver, A::SparseMatrixCSC{Float64, Int32}, b::Vector{Float64}, u::Vector{Float64}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::Bool, assump::OperatorAssumptions{Bool})
    @ LinearSolve ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:291
  [5] #init#3
    @ ~/.julia/packages/LinearSolve/O6GSm/src/common.jl:239 [inlined]
  [6] init
    @ ~/.julia/packages/LinearSolve/O6GSm/src/common.jl:140 [inlined]
  [7] #init#95
    @ ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:270 [inlined]
  [8] init
    @ ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:266 [inlined]
  [9] #construct_linear_solver#56
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/linear_solve.jl:77 [inlined]
 [10] construct_linear_solver
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/linear_solve.jl:61 [inlined]
 [11] init(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, alg::DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, J::SparseMatrixCSC{Float64, Int32}, fu::Vector{Float64}, u::Vector{Float64}; stats::SciMLBase.NLStats, pre_inverted::Val{false}, linsolve_kwargs::@NamedTuple{abstol::Float64, reltol::Float64}, abstol::Float64, reltol::Float64, timer::Nothing, alias_J::Bool, shared::Val{2}, kwargs::@Kwargs{})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/yZeYz/src/descent/damped_newton.jl:139
 [12] init
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/descent/damped_newton.jl:44 [inlined]
 [13] init(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, alg::GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, J::SparseMatrixCSC{Float64, Int32}, fu::Vector{Float64}, u::Vector{Float64}; shared::Val{1}, pre_inverted::Val{false}, linsolve_kwargs::@NamedTuple{abstol::Float64, reltol::Float64}, abstol::Float64, reltol::Float64, internalnorm::typeof(NonlinearSolveBase.L2_NORM), kwargs::@Kwargs{stats::SciMLBase.NLStats, timer::Nothing})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/yZeYz/src/descent/geodesic_acceleration.jl:83
 [14] __init(::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion{Float64}, GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, Nothing, Nothing, Nothing, Val{true}}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Nothing, reltol::Nothing, maxtime::Nothing, termination_condition::Nothing, internalnorm::Function, linsolve_kwargs::@NamedTuple{}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{})
    @ NonlinearSolveFirstOrder ~/.julia/packages/NonlinearSolveFirstOrder/UpLMt/src/solve.jl:167
 [15] __init
    @ ~/.julia/packages/NonlinearSolveFirstOrder/UpLMt/src/solve.jl:121 [inlined]
 [16] #__solve#146
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/solve.jl:5 [inlined]
 [17] __solve
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/solve.jl:1 [inlined]
 [18] #solve_call#36
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:667 [inlined]
 [19] solve_call
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:624 [inlined]
 [20] solve_up(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion{Float64}, GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, Nothing, Nothing, Nothing, Val{true}}; kwargs::@Kwargs{alias_u0::Bool})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1183
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1177 [inlined]
 [22] #solve#44
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1171 [inlined]
 [23] solve(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, args::GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion{Float64}, GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, Nothing, Nothing, Nothing, Val{true}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1136
 [24] top-level scope
    @ ~/Documents/julia/nonlinear-solve-mwe/mwe.jl:37

Environment (please complete the following information):

  • Output of using Pkg; Pkg.status()
Status `~/Documents/julia/nonlinear-solve-mwe/Project.toml`
  [8913a72c] NonlinearSolve v4.8.0
  [94fada2c] PowerFlows v0.9.0
  [f00506e0] PowerSystemCaseBuilder v1.3.11
  [bcd98974] PowerSystems v4.6.2
  [2f01184e] SparseArrays v1.11.0
  • Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
Status `~/Documents/julia/nonlinear-solve-mwe/Manifest.toml`
  [47edcb42] ADTypes v1.14.0
  [1520ce14] AbstractTrees v0.4.5
  [7d9f7c33] Accessors v0.1.42
  [79e6a3ab] Adapt v4.3.0
  [13e28ba4] AppleAccelerate v0.4.0
⌃ [4fba245c] ArrayInterface v7.18.0
  [4c555306] ArrayLayouts v1.11.1
  [62783981] BitTwiddlingConvenienceFunctions v0.1.6
  [a74b3585] Blosc v0.7.3
  [70df07ce] BracketingNonlinearSolve v1.2.0
  [2a0fbf3d] CPUSummary v0.2.6
  [336ed68f] CSV v0.10.15
  [d360d2e6] ChainRulesCore v1.25.1
  [fb6a15b2] CloseOpenIntervals v0.1.13
  [944b1d66] CodecZlib v0.7.8
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.1
  [f70d9fcc] CommonWorldInvalidations v1.0.0
  [34da2185] Compat v4.16.0
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.8
  [adafc99b] CpuId v0.3.1
  [a8cc5b0e] Crayons v4.1.1
  [a10d1c49] DBInterface v2.6.1
  [9a962f9c] DataAPI v1.16.0
  [a93c6f00] DataFrames v1.7.0
  [864edb3b] DataStructures v0.18.22
  [e2d170a0] DataValueInterfaces v1.0.0
  [8bb1440f] DelimitedFiles v1.9.1
⌃ [2b5f629d] DiffEqBase v6.170.1
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [a0c0ee7d] DifferentiationInterface v0.6.52
  [ffbed154] DocStringExtensions v0.9.4
  [4e289a0a] EnumX v1.0.5
  [f151be2c] EnzymeCore v0.8.8
  [e2ba6199] ExprTools v0.1.10
  [55351af7] ExproniconLite v0.10.14
  [7034ab61] FastBroadcast v0.3.5
  [9aa1b823] FastClosures v0.3.2
  [a4df4552] FastPower v1.1.2
  [48062228] FilePathsBase v0.9.24
  [1a297f60] FillArrays v1.13.0
  [6a86dc24] FiniteDiff v2.27.0
  [f6369f11] ForwardDiff v1.0.1
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [46192b85] GPUArraysCore v0.2.0
  [c8ec2601] H5Zblosc v0.1.2
  [f67ccb44] HDF5 v0.17.2
  [615f187c] IfElse v0.1.1
  [2cd47ed4] InfrastructureSystems v2.6.0
  [842dd82b] InlineStrings v1.4.3
  [3587e190] InverseFunctions v0.1.17
  [41ab1584] InvertedIndices v1.3.1
  [92d709cd] IrrationalConstants v0.2.4
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.7.0
  [0f8b85d8] JSON3 v1.14.2
  [ae98c720] Jieko v0.2.1
  [ef3ab10e] KLU v0.6.0
  [ba0b0d4f] Krylov v0.10.1
  [b964fa9f] LaTeXStrings v1.4.0
  [10f19ff3] LayoutPointers v0.1.17
  [5078a376] LazyArrays v2.6.1
  [1d6d02ad] LeftChildRightSiblingTrees v0.2.0
  [87fe0de2] LineSearch v0.1.4
⌃ [7ed4a6bd] LinearSolve v3.9.0
  [2ab3a3ac] LogExpFunctions v0.3.29
⌅ [33e6dc65] MKL v0.7.0
  [3da0fdf6] MPIPreferences v0.1.11
  [1914dd2f] MacroTools v0.5.16
  [d125e4d3] ManualMemory v0.1.8
  [bb5d69b7] MaybeInplace v0.1.4
  [e1d29d7a] Missings v1.2.0
  [2e0e35c7] Moshi v0.3.5
  [46d2c3a1] MuladdMacro v0.2.4
  [ffc61752] Mustache v1.0.20
  [77ba4419] NaNMath v1.1.3
  [8913a72c] NonlinearSolve v4.8.0
  [be0214bd] NonlinearSolveBase v1.6.0
  [5959db7a] NonlinearSolveFirstOrder v1.4.0
  [9a2c21bd] NonlinearSolveQuasiNewton v1.3.0
  [26075421] NonlinearSolveSpectralMethods v1.2.0
  [bac558e1] OrderedCollections v1.8.0
  [d96e819e] Parameters v0.12.3
⌅ [46dd5b70] Pardiso v0.5.7
  [69de0a69] Parsers v2.8.3
  [f517fe37] Polyester v0.7.16
  [1d0040c9] PolyesterWeave v0.2.2
  [2dfb63ee] PooledArrays v1.4.3
  [dd99e9e3] PowerFlowData v1.5.0
  [94fada2c] PowerFlows v0.9.0
⌅ [bed98974] PowerNetworkMatrices v0.12.1
  [f00506e0] PowerSystemCaseBuilder v1.3.11
  [bcd98974] PowerSystems v4.6.2
⌅ [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [08abe8d2] PrettyTables v2.4.0
  [33c8b6b6] ProgressLogging v0.1.4
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.33.0
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.14
  [94e857df] SIMDTypes v0.1.0
  [0aa819cd] SQLite v1.6.1
  [0bca4576] SciMLBase v2.86.2
  [19f34311] SciMLJacobianOperators v0.1.3
⌃ [c0aeaf25] SciMLOperators v0.3.13
  [53ae85a6] SciMLStructures v1.7.0
  [91c51154] SentinelArrays v1.4.8
  [efcf1570] Setfield v1.1.2
⌃ [727e6d20] SimpleNonlinearSolve v2.3.0
  [a2af1166] SortingAlgorithms v1.2.1
  [0a514795] SparseMatrixColorings v0.4.19
  [276daf66] SpecialFunctions v2.5.1
  [aedffcd0] Static v1.2.0
  [0d7ed370] StaticArrayInterface v1.8.0
  [1e83bf80] StaticArraysCore v1.4.3
  [10745b16] Statistics v1.11.1
  [7792a7ef] StrideArraysCore v0.5.7
  [69024149] StringEncodings v0.3.7
  [892a3eda] StringManipulation v0.4.1
  [856f2bd8] StructTypes v1.11.0
  [2efcf032] SymbolicIndexingInterface v0.3.40
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.0
  [5d786b92] TerminalLoggers v0.1.7
  [8290d209] ThreadingUtilities v0.5.3
  [9e3dc215] TimeSeries v0.24.2
  [a759f4b9] TimerOutputs v0.5.28
  [3bb67fe8] TranscodingStreams v0.11.3
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [ea10d353] WeakRefStrings v1.4.2
  [76eceee3] WorkerUtilities v1.6.1
  [ddb6d928] YAML v0.4.13
  [0b7ba130] Blosc_jll v1.21.7+0
  [0234f1f7] HDF5_jll v1.14.6+0
  [e33a78d0] Hwloc_jll v2.12.0+0
⌅ [1d5cc7b8] IntelOpenMP_jll v2024.2.1+0
  [94ce4f54] Libiconv_jll v1.18.0+0
  [5ced341a] Lz4_jll v1.10.1+0
⌅ [856f044c] MKL_jll v2024.2.0+0
  [7cb0a576] MPICH_jll v4.3.0+1
  [f1f71cc9] MPItrampoline_jll v5.5.3+0
  [9237b28f] MicrosoftMPI_jll v10.1.4+3
  [fe0851c0] OpenMPI_jll v5.0.7+2
  [458c3c95] OpenSSL_jll v3.5.0+0
  [efe28fd5] OpenSpecFun_jll v0.5.6+0
  [76ed43ae] SQLite_jll v3.48.0+0
  [3161d3a3] Zstd_jll v1.5.7+1
  [477f73a3] libaec_jll v1.1.3+0
  [1317d2d5] oneTBB_jll v2022.0.0+0
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [4af54fe1] LazyArtifacts v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.11.0
  [de0858da] Printf v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.6.0+0
  [e37daf67] LibGit2_jll v1.7.2+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.6+0
  [14a3606d] MozillaCACerts_jll v2023.12.12
  [4536629a] OpenBLAS_jll v0.3.27+1
  [05823500] OpenLibm_jll v0.8.5+0
  [bea87d4a] SuiteSparse_jll v7.7.0+0
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.59.0+0
  [3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
  • Output of versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 12 × Apple M2 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

Additional context

There's a chance that I'm just doing something wrong here, that this isn't actually a bug. But the error message is bizarre enough that I decided to pick "Bug" over "Question" when opening this issue. I also created a discourse post about this same issue.

Activity

luke-kiernan
added
bugSomething isn't working
on May 6, 2025
avik-pal

avik-pal commented on May 7, 2025

@avik-pal
avik-pal
on May 7, 2025
Member

LM needs to solve a linear system with a size 2N x N matrix. The bottom N x N part is a diagonal term. I am not sure why we are dispatching to KLU, so this might be a LinearSolve issue?

ChrisRackauckas

ChrisRackauckas commented on May 7, 2025

@ChrisRackauckas
ChrisRackauckas
on May 7, 2025
Member

https://github.com/SciML/LinearSolve.jl/blob/main/ext/LinearSolveSparseArraysExt.jl#L248

Can someone add an @show into there and see if assump.issq is true?

luke-kiernan

luke-kiernan commented on May 7, 2025

@luke-kiernan
luke-kiernan
on May 7, 2025 · edited by luke-kiernan
Author

Can someone add an @show into there and see if assump.issq is true?

I did so, and I'm getting that it's false. I also added a print statement to line 145 of the same file: indeed, the size of A is (20,10).

Edit: working my way up the stacktrace, the errant matrix is created in stack frame [12] or higher. i.e. if I add a print statement around line 138 of damped_newton.jl inside NonlinearSolveBase, the matrix A has wrong dimensions there too.

ChrisRackauckas

ChrisRackauckas commented on May 11, 2025

@ChrisRackauckas
ChrisRackauckas
on May 11, 2025
Member
using NonlinearSolve
using SciMLBase
using SparseArrays
using PowerFlows, PowerSystemCaseBuilder
const PF = PowerFlows
const PSB = PowerSystemCaseBuilder
# small system for demo purposes: actual target system has 10k buses. 
sys = PSB.build_system(PSB.MatpowerTestSystems, "matpower_case5_sys")
ERROR: KeyError: key :PowerSystems not found
Stacktrace:
  [1] root_module(where::Module, name::Symbol)
    @ Base ./loading.jl:1951
  [2] get_module(module_name::String)
    @ InfrastructureSystems ~/.julia/packages/InfrastructureSystems/yUXE4/src/utils/utils.jl:382
  [3] get_type_from_serialization_metadata(metadata::Dict{String, Any})
    @ InfrastructureSystems ~/.julia/packages/InfrastructureSystems/yUXE4/src/serialization.jl:132
  [4] get_type_from_serialization_data(data::Dict{String, Any})
    @ InfrastructureSystems ~/.julia/packages/InfrastructureSystems/yUXE4/src/serialization.jl:128
  [5] deserialize_components!(sys::PowerSystems.System, raw::Dict{String, Any})
    @ PowerSystems ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:1464
  [6] deserialize(::Type{…}, filename::String; time_series_read_only::Bool, time_series_directory::Nothing, kwargs::@Kwargs{…})
    @ PowerSystems ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:1441
  [7] deserialize
    @ ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:1389 [inlined]
  [8] PowerSystems.System(file_path::String; assign_new_uuids::Bool, kwargs::@Kwargs{})
    @ PowerSystems ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:204
  [9] _build_system(category::Type{…}, name::String, sys_descriptor::SystemDescriptor, print_stat::Bool; kwargs::@Kwargs{})
    @ PowerSystemCaseBuilder ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:112
 [10] _build_system
    @ ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:75 [inlined]
 [11] #build_system#236
    @ ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:26 [inlined]
 [12] build_system
    @ ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:18 [inlined]
 [13] build_system(category::Type{MatpowerTestSystems}, name::String)
    @ PowerSystemCaseBuilder ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:18
 [14] top-level scope
    @ ~/Desktop/test2.jl:8
Some type information was truncated. Use `show(err)` to see complete types.
ChrisRackauckas

ChrisRackauckas commented on May 12, 2025

@ChrisRackauckas
ChrisRackauckas
on May 12, 2025
Member

And I couldn't recreate it with a simple MWE:

using NonlinearSolve

function f!(du, u, p)
    du[1] = 2u[1] - 2
    du[2] = u[1] - 4u[2]
end

jac = SparseMatrixCSC([
    1 0
    1 1
])

u0 = [0.0, 0.0]
prob = NonlinearProblem(
    NonlinearFunction(f!, jac_prototype = jac), u0)
solver = NonlinearSolveFirstOrder.LevenbergMarquardt()
solve(prob, solver)

Please try and make an MWE that doesn't require the extra packages. The core issue shouldn't require anything specific about f, it's just how you setup the problem, so it would be good to just choose a simple f and recreate the rest.

luke-kiernan

luke-kiernan commented on May 12, 2025

@luke-kiernan
luke-kiernan
on May 12, 2025
Author

ERROR: KeyError: key :PowerSystems not found

Try adding PowerSystems as well.

And I couldn't recreate it with a simple MWE:

Strangely, if I change the indexing type to Int32, then I get the same error:

using NonlinearSolve
using SparseArrays

function f!(du, u, p)
    du[1] = 2u[1] - 2
    du[2] = u[1] - 4u[2]
end

jac = SparseMatrixCSC{Float64, Int32}([
    1.0 0.0
    1.0 1.0
])

u0 = [0.0, 0.0]
prob = NonlinearProblem(
    NonlinearFunction(f!, jac_prototype = jac), u0)
solver = NonlinearSolveFirstOrder.LevenbergMarquardt()
solve(prob, solver)
ChrisRackauckas
added a commit that references this issue on May 12, 2025

Fix SparseMatrixCSC with alternative Integer types

e2b9740
ChrisRackauckas
mentioned this on May 12, 2025
ChrisRackauckas
added a commit that references this issue on May 12, 2025

Fix non-square sparse matrix preallocations for defaults

d763a37
ChrisRackauckas
mentioned this on May 12, 2025
ChrisRackauckas
closed this as completedin #462on May 12, 2025
ChrisRackauckas

ChrisRackauckas commented on May 12, 2025

@ChrisRackauckas
ChrisRackauckas
on May 12, 2025
Member

Thanks, that pinpointed where the issue was. ArrayInterface and LinearSolve need to handle Int32 cases, so it'll need to wait for the LinearSolve.jl bump but this is fixed locally and just needs one more merge.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

      Development

      Participants

      @ChrisRackauckas@avik-pal@luke-kiernan

      Issue actions
        Levenberg-Marquardt with Sparse J: J isn't square?! · Issue #599 · SciML/NonlinearSolve.jl