# Maximum of matrix vector product for "large" matrices

**URL:** <https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229>\
**Category:** General Usage\
**Tags:** modelingtoolkit\
**Created:** [October 25, 2022, 3:35am UTC](https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229 "2022-10-25T03:35:22Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![DorianC](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dorianc/32/33792_2.png) [@DorianC](https://discourse.julialang.org/u/DorianC)\
**Post date:** [October 25, 2022, 3:35am UTC](https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229/1 "2022-10-25T03:35:22Z")

</div>

I am using MTK and noticed that when I construct an algebraic variable defined as the max of the matrix-vector product of the state vector, the results are accurate for small matrices, but go inaccurate rapidly as the matrix size increases.  
Specifically, I get accurate results up 100x10. For 500x10, it’s always inaccurate.  
Copy-paste example (you can simply uncomment the definition of C):

```julia
using LinearAlgebra
using ModelingToolkit, Plots, DifferentialEquations
using Symbolics: scalarize

@variables t
D = Differential(t)

function F(;name, α0 = fill(0.0, 10), C = fill(1.0, 10, 10))

    sts = @variables T(t)=0.0 α(t)[1:length(α0)]=α0

    eqs = [ scalarize(D.(α) .~ 0.0)
            T ~ maximum(C*α)
            ]

    ODESystem(eqs, t, [α..., T], []; name=name)
end

α0 = rand(10)

#C = rand(100, 10)
C = rand(500, 10)

@named mySys = F(α0=α0, C=C)

sys = structural_simplify(mySys)

prob = ODEProblem(sys, [], (0.0, 1.0))
sol = solve(prob)

println(maximum(C*α0))
println(sol[mySys.T][end])

```

Any clues how to get around that using MTK? If i build the same system directly in DifferentialEquations with an algebraic state defined the same manner, I don’t run into the same problem.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 25, 2022, 3:43am UTC](https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229/2 "2022-10-25T03:43:12Z")

</div>

@YingboMa @shashi could this be fixed by using the MuladdMacro expression changing on the generated code?

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [October 28, 2022, 7:59pm UTC](https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229/3 "2022-10-28T19:59:19Z")

</div>

I see very similar values for the last two lines

```julia
julia> println(maximum(C*α0))
4.421761822317552

julia> println(sol[mySys.T][end])
4.421761822317553

```

I started a bigger system (1000x100) running but it hasn’t finished yet…

```julia
(MTK) pkg> st
Project MTK v0.1.0
Status `~/Desktop/Julia/MTK/Project.toml`
  [0c46a032] DifferentialEquations v7.6.0
  [961ee093] ModelingToolkit v8.30.0
  [91a5bcdd] Plots v1.35.5
  [731186ca] RecursiveArrayTools v2.32.1
  [0c5d862f] Symbolics v4.13.0

```

Edit  
1000x100 finished!

```julia
julia> println(maximum(C*α0))
33.386913594597175

julia> println(sol[mySys.T][end])
33.38691359459718

```

---

<div class="post-metadata">

**Author:** ![DorianC](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dorianc/32/33792_2.png) [@DorianC](https://discourse.julialang.org/u/DorianC)\
**Post date:** [October 30, 2022, 8:21pm UTC](https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229/4 "2022-10-30T20:21:18Z")

</div>

Ah, good find! I checked and made sure the package versions were the same, but it didn’t make a difference in Julia 1.7.2. When I tried in Julia 1.8, I got nearly the same result (not equal but good enough)! Thank you for pointing that out!

---

<div class="post-metadata">

**Author:** ![DorianC](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dorianc/32/33792_2.png) [@DorianC](https://discourse.julialang.org/u/DorianC)\
**Post date:** [October 31, 2022, 8:04pm UTC](https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229/5 "2022-10-31T20:04:53Z")

</div>

Ok, now for part II now that the max issue is resolved by switching to Julia 1.8.2 🙂 :

My original code looks a lot more like:

```julia
using LinearAlgebra
using ModelingToolkit, Plots, DifferentialEquations
using Symbolics: scalarize

@variables t
D = Differential(t)

function myComponent(;name, Nm = 10, α0=zeros(Nm))

    sts = @variables α(t)[1:Nm]=α0
    eqs = scalarize(D.(α) .~ 0.0)

    ODESystem(eqs, t, [α...], []; name)
end

function F(;name, α0s = [fill(0.0, 10) for i in 1:20], C = fill(1.0, 10, 200))

    components = [myComponent(name = Symbol(string(i)), α0 = α0s[i] ) for i in 1:length(α0s)]

    α = reduce(vcat, [component.α for component in components])

    @variables T(t)=0.0

    eqs = [ T ~ maximum(C*α)
            ]

    compose(ODESystem(eqs, t, [T], [],; name), components...)

end

α0 = fill(rand(10), 20)

#C = rand(100, 10)
C = rand(500, 200)

@named mySys = F(α0s=α0, C=C)

sys = structural_simplify(mySys)

prob = ODEProblem(sys, [], (0.0, 1.0))
@time solve(prob)

```

When I run that on Julia 1.7.2. I get:

```julia
5.666656 seconds (10.51 M allocations: 634.450 MiB, 1.81% gc time, 99.94% compilation time)
retcode: Success

```

When i try to run that on Julia 1.8.2, it takes \> 5min (I’m still waiting for the answer and will edit).

Now I’m sure I’m doing something dumb with this concatenation (if this wasn’t symbolic, i’d at least pre-allocate). Any tips for a newbie at symbolic math/MTK?

Thanks for the help!

---

<div class="post-metadata">

**Author:** ![dorian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dorian/32/43827_2.png) [@dorian](https://discourse.julialang.org/u/dorian)\
**Post date:** [October 31, 2022, 8:17pm UTC](https://discourse.julialang.org/t/maximum-of-matrix-vector-product-for-large-matrices/89229/6 "2022-10-31T20:17:06Z")

</div>

> [@DorianC](#):
>
> it takes \> 5min (I’m still waiting for the answer and will edit).

It’s been 15min now, so I give up.