# Speed up solution of large system of equations
**URL:** https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659
**Category:** Performance
**Tags:** nlsolve, nonlinear, speed-optimization
**Created:** [August 10, 2023, 2:25am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659 "2023-08-10T02:25:12Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 2:25am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/1 "2023-08-10T02:25:12Z")
I have a system of equations:
\begin{aligned} & S^m\_{a\_m, e\_m}-\mu^{0m}\_{a\_m, e\_m}=\sum\_{a\_m=1}^{46} {\Pi}\_{a\_m, e\_m, a\_f, e\_f}\sqrt{S^m\_{a\_m, e\_m} S^f\_{a\_f, e\_f}} \prod\_{k=0}^{T(a\_m, a\_f)-1}\left(\frac{\mu^{0m}\_{a\_m+k, e\_m} \mu^{0f}\_{a\_f+k, e\_f}}{S^m\_{a\_m+k, e\_m} S^f\_{a\_f+k, e\_f}}\right)^{(1 / 2)(\beta(1-\delta))^k}, \\ & S^f\_{a\_f, e\_f}-\mu^{0f}\_{a\_f, e\_f}=\sum\_{a\_f=1}^{46} {\Pi}\_{a\_m, e\_m, a\_f, e\_f}\sqrt{S^m\_{a\_m, e\_m} S^f\_{a\_f, e\_f}} \prod\_{k=0}^{T(a\_m, a\_f)-1}\left(\frac{\mu^{0m}\_{a\_m+k, e\_m} \mu^{0f}\_{a\_f+k, e\_f}}{S^m\_{a\_m+k, e\_m} S^f\_{a\_f+k, e\_f}}\right)^{(1 / 2)(\beta(1-\delta))^k}. \end{aligned}
for a\_m, a\_f \in \{1, \ldots, 46\} and e\_m, e\_f \in \{1, 2\}, where T(a\_m, a\_f) = 46 - \max\{a\_m, a\_f\} + 1. S^m, S^f, \Pi, \mu^{0m}, \mu^{0f} are matrices. I need to solve that system for each of \mu^{0m}\_{a\_m, e\_m}, \mu^{0f}\_{a\_f, e\_f}.
I have a script that solves that system using a generalization of the [IPFP](https://en.wikipedia.org/wiki/Iterative_proportional_fitting) algorithm:
```julia
using NLsolve
using LinearAlgebra
const Z = 46
const δ = 0.1
const β = 0.95
T(am, af) = Z - max(am, af) + 1
const IDX_ae = CartesianIndices((1:46, 1:2))
function insideproduct_num(am, em, af, ef; μ0m=μ0m, μ0f=μ0f)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(μ0m), eltype(μ0f)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(abs.(μ0m[am + k, em] * μ0f[af + k, ef]))
power *= exponent
end
return exp(p)
end
function insideproduct_den(am, em, af, ef; Sm=Sm, Sf=Sf)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(Sm), eltype(Sf)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(abs.(Sm[am + k, em] * Sf[af + k, ef]))
power *= exponent
end
return exp(p)
end
function fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f)
for idx in CartesianIndices(auxmat)
am, em, af, ef = Tuple(idx)
auxmat[am, em, af, ef] = Π[am, em, af, ef] * sqrt(Sm[am, em]*Sf[af, ef]) /
insideproduct_den(am, em, af, ef; Sm=Sm, Sf=Sf)
end
return auxmat
end
@views function MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f)
for idx_ae in IDX_ae
am, em = Tuple(idx_ae)
s = zero(eltype(μ0m))
for idx_fi in IDX_ae
s += auxmat[am, em, idx_fi[1], idx_fi[2]] * insideproduct_num(am, em, idx_fi[1], idx_fi[2]; μ0m=μ0m, μ0f=μ0f)
end
res[am, em] = Sm[am, em] - μ0m[am, em] - s
end
return res
end
@views function MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m)
for idx_ae in IDX_ae
af, ef = Tuple(idx_ae)
s = zero(eltype(μ0f))
for idx_mi in IDX_ae
s += auxmat[idx_mi[1], idx_mi[2], af, ef] * insideproduct_num(idx_mi[1], idx_mi[2], af, ef; μ0m=μ0m, μ0f=μ0f)
end
res[af, ef] = Sf[af, ef] - μ0f[af, ef] - s
end
return res
end
function ipfp(Π, Sm, Sf; initμ0=initμ0, auxmat=auxmat, tol=1, method=:trust_region, autodiff=:forward, iterations=1_000,
innerxtol=1, innerftol=1)
i = 0
μ0f1 = Sf
μ0m = initμ0
d = Inf
while d > tol && i <= iterations
fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f1)
solM = NLsolve.nlsolve((res, μ0m) -> MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f1), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0m = solM.zero
solF = NLsolve.nlsolve((res, μ0f) -> MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0f2 = solF.zero
d = norm(μ0f1 - μ0f2, Inf)
μ0f1 = μ0f2
i += 1
end
println("i = ", i, ", d = ", d)
return abs.(μ0m), abs.(μ0f1)
end
```
Initialize with
```julia
Π = 5 .* rand(46, 2, 46, 2)
Sm = 5 .* rand(46, 2)
Sf = 5 .* rand(46, 2)
μ0m = similar(Sm)
μ0f = similar(Sf)
auxmat = similar(Π)
initμ0 = 5 .* ones(size(Sm))
```
Now run it
```julia
ipfp(Π, Sm, Sf; initμ0= initμ0)
```
I would like to speed that up as much as possible. Right now it takes over 2 seconds to solve this toy example. Are there any obvious speed ups I’m missing? Perhaps multithreading?
Note: If I scale the matrices by a factor larger than 5, the algorithm slows down even more.
---
### Author:  [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)
#### Post date: [August 10, 2023, 4:54am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/2 "2023-08-10T04:54:50Z")
So in situations like this I always check three things to see what might be slowing the code down:
- `@code_warntype` to see if there are problems with identifying variable types
- `@benchmark` from the BenchmarkTools.jl package to look at allocations
- `@profview` to see where the function is spending most of its time.
When I ran `@code_warntype ipfp(Π, Sm, Sf; initμ0= initμ0)` I noticed that `auxmat` came back as type Any.
I think the issue is you have to pass the variables down through each function. When you write:
```julia
function ipfp(Π, Sm, Sf; initμ0=initμ0, auxmat=auxmat) #plus other keywords
...
end
```
The function `ipfp` doesn’t know what auxmat is because it hasn’t been defined yet. So when you call it without giving it auxmat, it has to go looking for the variable in the global scope.
You could instead write:
```julia
function ipfp(Π, Sm, Sf; initμ0=initμ0, auxmat) #no default value for auxmat
...
end
```
and call the function as:
```julia
ipfp(Π, Sm, Sf; initμ0=initμ0, auxmat=auxmat)
#This works too
#Julia is cool and will match keywords with variable names
ipfp(Π, Sm, Sf; initμ0, auxmat)
```
This will have to be fixed for the other function definitions as well.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 5:40am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/3 "2023-08-10T05:40:07Z")
Thank you for that. I wasn’t aware of that problem. I still get some Any types. I suspect the problem is that Julia doesn’t know the type of the object that returns the zero of the system inside `ipfp`. Any suggestions for that?
```julia
using NLsolve
using LinearAlgebra
using Random
Random.seed!(1293804)
const Z = 46
const δ = 0.1
const β = 0.95
T(am, af) = Z - max(am, af) + 1
const IDX_ae = CartesianIndices((1:46, 1:2))
function insideproduct_num(am, em, af, ef; μ0m, μ0f)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(μ0m), eltype(μ0f)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(abs.(μ0m[am + k, em] * μ0f[af + k, ef]))
power *= exponent
end
return exp(p)
end
function insideproduct_den(am, em, af, ef; Sm, Sf)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(Sm), eltype(Sf)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(abs.(Sm[am + k, em] * Sf[af + k, ef]))
power *= exponent
end
return exp(p)
end
function fillauxmat!(auxmat; Π, Sm, Sf)
for idx in CartesianIndices(auxmat)
am, em, af, ef = Tuple(idx)
auxmat[am, em, af, ef] = Π[am, em, af, ef] * sqrt(Sm[am, em]*Sf[af, ef]) /
insideproduct_den(am, em, af, ef; Sm=Sm, Sf=Sf)
end
return auxmat
end
@views function MMsysM!(res, μ0m; auxmat, Sm, μ0f)
for idx_ae in IDX_ae
am, em = Tuple(idx_ae)
s = zero(eltype(μ0m))
for idx_fi in IDX_ae
s += auxmat[am, em, idx_fi[1], idx_fi[2]] * insideproduct_num(am, em, idx_fi[1], idx_fi[2]; μ0m=μ0m, μ0f=μ0f)
end
res[am, em] = Sm[am, em] - μ0m[am, em] - s
end
return res
end
@views function MMsysF!(res, μ0f; auxmat, Sf, μ0m)
for idx_ae in IDX_ae
af, ef = Tuple(idx_ae)
s = zero(eltype(μ0f))
for idx_mi in IDX_ae
s += auxmat[idx_mi[1], idx_mi[2], af, ef] * insideproduct_num(idx_mi[1], idx_mi[2], af, ef; μ0m=μ0m, μ0f=μ0f)
end
res[af, ef] = Sf[af, ef] - μ0f[af, ef] - s
end
return res
end
function ipfp(Π, Sm, Sf; initμ0, auxmat, tol=1e-3, method=:trust_region, autodiff=:forward, iterations=1_000,
innerxtol=1e-3, innerftol=1e-3)
i = 0
μ0f1 = Sf
μ0m = initμ0
d = Inf
fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf)
while d > tol && i <= iterations
solM = NLsolve.nlsolve((res, μ0m) -> MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f1), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0m = solM.zero
solF = NLsolve.nlsolve((res, μ0f) -> MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0f2 = solF.zero
d = norm(μ0f1 - μ0f2, Inf)
μ0f1 = μ0f2
i += 1
end
println("i = ", i, ", d = ", d)
return abs.(μ0m), abs.(μ0f1)
end
```
---
### Author:  [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)
#### Post date: [August 10, 2023, 6:39am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/4 "2023-08-10T06:39:34Z")
Here is a slightly better version that doesn’t given Any types for the output:
```julia
function ipfp(Π, Sm, Sf; initμ0, auxmat, tol=1, method=:trust_region, autodiff=:forward, iterations=1_000, innerxtol=1, innerftol=1)
i = 0
μ0f1 = Sf #maybe μ0f1 = copy(Sf) instead?
μ0f2 = Sf
μ0m = initμ0
d = Inf
while d > tol && i <= iterations
fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f1)
solM = NLsolve.nlsolve((res, μ0m) -> MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f1), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0m .= solM.zero
solF = NLsolve.nlsolve((res, μ0f) -> MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0f2 .= solF.zero
d = norm(μ0f1 - μ0f2, Inf)
μ0f1 .= μ0f2
i += 1
end
println("i = ", i, ", d = ", d)
return abs.(μ0m), abs.(μ0f1)
end
```
I think there is some confusion when you assign one variable to another. When you write `μ0f1 = Sf` you have two variables pointing to the same data in memory (updating `μ0f1` will also update `Sf`). If you don’t want Sf to change you can write `μ0f1 = copy(Sf)`. Also here I use the `.=` notation which takes all the values on the right and assigns them to the same positions on the left (assuming the right and left variables are already defined).
There’s still something weird happening with the `@code_warntype` output. It looks like a temporary variable is being defined or something?
```julia
@code_warntype ipfp(Π, Sm, Sf; initμ0, auxmat)
...
Locals
@_7::Any
...
```
I think it has something to do with the nlsolve call but I’m not sure.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 10, 2023, 9:04am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/5 "2023-08-10T09:04:54Z")
If you profile, what’s taking all of the time? I presume setting up multiple nonlinear solves is not the way to go. Also, generally for large problems using NonlinearSolve and specializing the linear solver towards something for large (sparse) Jacobians is required for anything large.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 10:29am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/6 "2023-08-10T10:29:27Z")
`@profview` doesn’t seem to be working, the pane shows up blank. I have tried using NonlinearSolve, but I can’t figure out how to specify the problem ☹
---
### Author:  [@ctkelley](https://discourse.julialang.org/u/ctkelley)
#### Post date: [August 10, 2023, 11:15am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/7 "2023-08-10T11:15:08Z")
While it would be some extra work, Newton-Krylov might help you. There’s some code for that [here](https://github.com/ctkelley/SIAMFANLEquations.jl).
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 10, 2023, 11:24am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/8 "2023-08-10T11:24:50Z")
> [@amrods](#):
>
> I have tried using NonlinearSolve, but I can’t figure out how to specify the problem ☹
What did you try? It’s almost exactly the same as NLsolve in the high level interface, you just add `p` (which you can then reuse to remake).
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [August 10, 2023, 5:15pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/9 "2023-08-10T17:15:06Z")
> [@amrods](#):
>
> `@profview` doesn’t seem to be working, the pane shows up blank
This used to be a [bug in the VSCode extension](https://github.com/julia-vscode/julia-vscode/issues/3316) but it was fixed recently. Upgrading the extension should suffice, but if not, does the flame graph appear when you re-run the profiling command a second time?
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 6:55pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/10 "2023-08-10T18:55:51Z")
I wanted to reuse the work I have done. In particular the functions `MMsysM!` and `MMsysF!`.
Those functions are called as `MMsysM!(res, μ0m; auxmat, Sm, μ0f)`, where the kw arguments are basically parameters, however `auxmat` and `Sm` don’t change for a particular problem.
I did this:
```julia
probM = NonlinearProblem{isinplace}((res, μ0m, μ0f) -> MMsysM!(res, μ0m; auxmat, Sm, μ0f), initμ0, μ0f)
```
I am not sure whether the mutating argument is supposed to be the first one. Also, I am not sure how to change the parameters after `probM` is defined. I couldn’t find that in the documentation.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 7:00pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/11 "2023-08-10T19:00:04Z")
I didn’t know I had to run it twice. I was giving up too easily. How can I export the graph so I can upload it here? (And how do I read interpret it). [Here](https://github.com/amrods/sandbox/blob/main/profile.html) is a link to the html file.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [August 10, 2023, 7:30pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/12 "2023-08-10T19:30:22Z")
> [@amrods](#):
>
> I didn’t know I had to run it twice.
You shouldn’t have to run it twice, that’s why it’s a bug 😉
> [@amrods](#):
>
> How can I export the graph so I can upload it here?
I’m not sure how you did it, but it looks like you managed just fine by linking to the HTML. It’s much more useful than a static screenshot anyway because it’s interactive.
And in VSCode it’s even nicer because when you control + click the profile it leads you to the relevant line in your code.
[https://www.julia-vscode.org/docs/stable/userguide/profiler/](https://www.julia-vscode.org/docs/stable/userguide/profiler/)
> [@amrods](#):
>
> And how do I read interpret it
Roughly speaking, you want the tiles to be:
- mostly blue, because yellow = allocations (your code uses a lot of memory) and red = runtime dispatch (your code struggles to infer types)
- related to the functions where you expect to be spending time (that’s problem-specific)
At first glance, your profile doesn’t scream “underoptimized” to me, because most of the time is spent in mathematical operations, and the red tiles are related to autodiff (thus sometimes unavoidable). But I don’t know the details.
To go further, check out this blog post, which is the best resource I know on optimizing Julia code. It has a profiling section and much more:
[https://viralinstruction.com/posts/optimise/#how\_to\_optimise\_julia\_code\_a\_practical\_guide](https://viralinstruction.com/posts/optimise/#how_to_optimise_julia_code_a_practical_guide)
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 7:36pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/13 "2023-08-10T19:36:04Z")
It seems the fix for the bug hasn’t found its way into a release. I have version 1.47.2. I didn’t want to install a pre-release version. Thanks for the resources.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 10, 2023, 8:02pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/14 "2023-08-10T20:02:47Z")
> [@amrods](#):
>
> I am not sure whether the mutating argument is supposed to be the first one. Also, I am not sure how to change the parameters after `probM` is defined. I couldn’t find that in the documentation.
It looks like just `μ0f` is your parameters? Then you can just enclose the others
```julia
probM = NonlinearProblem{isinplace}((res, μ0f) -> MMsysM!(res, μ0m; auxmat, Sm, μ0f), initμ0, μ0f)
```
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 9:50pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/15 "2023-08-10T21:50:54Z")
I don’t need to specify the arguments? In the example in the documentation f(u, p), u is the input of the function, and p the parameters. That gets translated to `probN = NonlinearProblem(f, u0, p)`.
In my case, `μ0m` is the input, `μ0f` the parameters, but `res` is the array where the residuals of the system are stored. Just to clarify, writing `isinplace` takes care of all that so that I can specify the function as you did?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 10, 2023, 10:05pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/16 "2023-08-10T22:05:37Z")
Do you have an MWE? I’m trying to run the following:
```julia
using NLsolve
using LinearAlgebra
using Random
Random.seed!(1293804)
const Z = 46
const δ = 0.1
const β = 0.95
T(am, af) = Z - max(am, af) + 1
const IDX_ae = CartesianIndices((1:46, 1:2))
function insideproduct_num(am, em, af, ef; μ0m, μ0f)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(μ0m), eltype(μ0f)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(abs.(μ0m[am + k, em] * μ0f[af + k, ef]))
power *= exponent
end
return exp(p)
end
function insideproduct_den(am, em, af, ef; Sm, Sf)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(Sm), eltype(Sf)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(abs.(Sm[am + k, em] * Sf[af + k, ef]))
power *= exponent
end
return exp(p)
end
function fillauxmat!(auxmat; Π, Sm, Sf)
for idx in CartesianIndices(auxmat)
am, em, af, ef = Tuple(idx)
auxmat[am, em, af, ef] = Π[am, em, af, ef] * sqrt(Sm[am, em]*Sf[af, ef]) /
insideproduct_den(am, em, af, ef; Sm=Sm, Sf=Sf)
end
return auxmat
end
@views function MMsysM!(res, μ0m; auxmat, Sm, μ0f)
for idx_ae in IDX_ae
am, em = Tuple(idx_ae)
s = zero(eltype(μ0m))
for idx_fi in IDX_ae
s += auxmat[am, em, idx_fi[1], idx_fi[2]] * insideproduct_num(am, em, idx_fi[1], idx_fi[2]; μ0m=μ0m, μ0f=μ0f)
end
res[am, em] = Sm[am, em] - μ0m[am, em] - s
end
return res
end
@views function MMsysF!(res, μ0f; auxmat, Sf, μ0m)
for idx_ae in IDX_ae
af, ef = Tuple(idx_ae)
s = zero(eltype(μ0f))
for idx_mi in IDX_ae
s += auxmat[idx_mi[1], idx_mi[2], af, ef] * insideproduct_num(idx_mi[1], idx_mi[2], af, ef; μ0m=μ0m, μ0f=μ0f)
end
res[af, ef] = Sf[af, ef] - μ0f[af, ef] - s
end
return res
end
function ipfp(Π, Sm, Sf; initμ0, auxmat, tol=1e-3, method=:trust_region, autodiff=:forward, iterations=1_000,
innerxtol=1e-3, innerftol=1e-3)
i = 0
μ0f1 = Sf
μ0m = initμ0
d = Inf
fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf)
while d > tol && i <= iterations
solM = NLsolve.nlsolve((res, μ0m) -> MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f1), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0m = solM.zero
solF = NLsolve.nlsolve((res, μ0f) -> MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0f2 = solF.zero
d = norm(μ0f1 - μ0f2, Inf)
μ0f1 = μ0f2
i += 1
end
println("i = ", i, ", d = ", d)
return abs.(μ0m), abs.(μ0f1)
end
Π = 5 .* rand(46, 2, 46, 2)
Sm = 5 .* rand(46, 2)
Sf = 5 .* rand(46, 2)
μ0m = similar(Sm)
μ0f = similar(Sf)
auxmat = similar(Π)
initμ0 = 5 .* ones(size(Sm))
ipfp(Π, Sm, Sf; initμ0, auxmat)
```
but it just runs forever.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 10, 2023, 10:44pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/17 "2023-08-10T22:44:38Z")
Sorry about that. This should work:
```julia
using NLsolve
using LinearAlgebra
using Random
Random.seed!(1293804)
const Z = 46
const δ = 0.1
const β = 0.95
T(am, af) = Z - max(am, af) + 1
const IDX_ae = CartesianIndices((1:46, 1:2))
function insideproduct_num(am, em, af, ef; μ0m, μ0f)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(μ0m), eltype(μ0f)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(abs.(μ0m[am + k, em] * μ0f[af + k, ef]))
power *= exponent
end
return exp(p)
end
function insideproduct_den(am, em, af, ef; Sm, Sf)
exponent = β*(1 - δ)
power = 1.0
p = zero(promote_type(eltype(Sm), eltype(Sf)))
for k in 0:T(am, af)-1
p += 0.5 * power * log(Sm[am + k, em] * Sf[af + k, ef])
power *= exponent
end
return exp(p)
end
function fillauxmat!(auxmat; Π=Π, Sm, Sf)
for idx in CartesianIndices(auxmat)
am, em, af, ef = Tuple(idx)
auxmat[am, em, af, ef] = Π[am, em, af, ef] * sqrt(Sm[am, em]*Sf[af, ef]) /
insideproduct_den(am, em, af, ef; Sm=Sm, Sf=Sf)
end
return auxmat
end
function MMsysM!(res, μ0m; auxmat, Sm, μ0f)
for idx_ae in IDX_ae
am, em = Tuple(idx_ae)
s = zero(eltype(μ0m))
for idx_fi in IDX_ae
s += auxmat[am, em, idx_fi[1], idx_fi[2]] * insideproduct_num(am, em, idx_fi[1], idx_fi[2]; μ0m=μ0m, μ0f=μ0f)
end
res[am, em] = Sm[am, em] - μ0m[am, em] - s
end
return res
end
function MMsysF!(res, μ0f; auxmat, Sf, μ0m)
for idx_ae in IDX_ae
af, ef = Tuple(idx_ae)
s = zero(eltype(μ0f))
for idx_mi in IDX_ae
s += auxmat[idx_mi[1], idx_mi[2], af, ef] * insideproduct_num(idx_mi[1], idx_mi[2], af, ef; μ0m=μ0m, μ0f=μ0f)
end
res[af, ef] = Sf[af, ef] - μ0f[af, ef] - s
end
return res
end
function ipfp(Π, Sm, Sf; initμ0, auxmat, tol=1, method=:trust_region, autodiff=:forward, iterations=1_000,
innerxtol=1, innerftol=1)
i = 0
μ0f1 = copy(Sf)
μ0f2 = copy(Sf)
μ0m = initμ0
d = Inf
fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf)
while d > tol && i <= iterations
solM = NLsolve.nlsolve((res, μ0m) -> MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f1), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0m .= solM.zero
solF = NLsolve.nlsolve((res, μ0f) -> MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m), initμ0; method=method, autodiff=autodiff,
xtol=innerxtol, ftol=innerftol, iterations=iterations)
μ0f2 .= solF.zero
d = norm(μ0f1 - μ0f2, Inf)
μ0f1 .= μ0f2
i += 1
end
print("iter = ", i, ", tol = ", d)
return abs.(μ0m), abs.(μ0f1)
end
Π1 = rand(46, 2, 46, 2)
Sm1 = rand(46, 2)
Sf1 = rand(46, 2)
μ0m = similar(Sm1)
μ0f = similar(Sf1)
auxmat = zeros(size(Π1))
initμ0 = 1 .* ones(size(Sm1))
res1 = similar(initμ0)
@time sol1 = ipfp(Π1, Sm1, Sf1; initμ0=initμ0, auxmat=auxmat, tol=1e-4)
```
I’m trying to test the somewhat full system, but with toy inputs to assess speed.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 13, 2023, 1:16am UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/18 "2023-08-13T01:16:04Z")
Two further questions:
1. I noticed that solving the system slows down for “larger” values of S\_m and S\_f. For example:
```julia
Sm5 = 5 .* Sm1
Sf5 = 5 .* Sf1
```
Since the system is homogeneous of degree zero in S\_m, S\_f, \mu^{0m}, \mu^{0f}, is it ok to scale down the system, solve it (presumably quicker) and then scale it back up? I note that all those variables must be in the same units.
1. The way this generalization of IPFP works is by sequentially solving subsystems until convergence. I thought of having the tolerance of solving those subsystems depend on how close the global algorithm is from convergence. That way, when the global algorithm is far from convergence, the subsystems do not need to be solved to high precision (and so will solve quickly), but when the global algorithm is close to convergence, then subsystems will be solved to higher precision. I modified the `ipfp` function to:
```julia
function ipfp_ad(Π, Sm, Sf; initμ0, auxmat, tol=1, method=:trust_region, autodiff=:forward, iterations=1_000)
i = 0
μ0f1 = copy(Sf)
μ0f2 = copy(Sf)
μ0m = initμ0
d = Inf
fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf)
while d > tol && i <= iterations
if d > 1
solM = NLsolve.nlsolve((res, μ0m) -> MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f1), initμ0; method=method, autodiff=autodiff,
xtol=1, ftol=1, iterations=iterations)
μ0m .= solM.zero
solF = NLsolve.nlsolve((res, μ0f) -> MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m), initμ0; method=method, autodiff=autodiff,
xtol=1, ftol=1, iterations=iterations)
μ0f2 .= solF.zero
d = norm(μ0f1 - μ0f2, Inf)
μ0f1 .= μ0f2
else
solM = NLsolve.nlsolve((res, μ0m) -> MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f1), initμ0; method=method, autodiff=autodiff,
xtol=d, ftol=d, iterations=iterations)
μ0m .= solM.zero
solF = NLsolve.nlsolve((res, μ0f) -> MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m), initμ0; method=method, autodiff=autodiff,
xtol=d, ftol=d, iterations=iterations)
μ0f2 .= solF.zero
d = norm(μ0f1 - μ0f2, Inf)
μ0f1 .= μ0f2
end
i += 1
end
print("iter = ", i, ", tol = ", d)
return abs.(μ0m), abs.(μ0f1)
end
```
What do people think of those 2 things?
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [August 19, 2023, 3:19pm UTC](https://discourse.julialang.org/t/speed-up-solution-of-large-system-of-equations/102659/19 "2023-08-19T15:19:48Z")
anyone?