Question about determining solution set matrices (D_v, D_i, ...) regarding exact arithmetic

[hemmer]

Thank you very much for making this library available, it's been very instructive for me learning about state-space methods. As a learning excerise, I've been re-implementing myself using python/sympy (a symbolic analysis library). I've started by implemented the linear version (just for some simple RC filters) using discrete time system eqn 13a in [1]. This is straightforward as there is a unique solution for the linear matrix problem, and sympy has a symbolic LUsolve routine so I can solve this symbolically (i.e. without component values being subbed in yet). This works as expected, and I've verified against analytical solutions, other solvers etc.

Now thinking about the nonlinear case, I've been reading the code for gensolve routine. If I understand, BigInt is used to avoid worrying about round-off errors during the initial solve (i.e. when gensolve is called in the setup phase from model_matrices(...)). However in the gensolve routine, there is still apparently a threshold/tolerance generated, and a reliance on floating point arithmetic.

https://github.com/HSU-ANT/ACME.jl/blob/main/src/ACME.jl#L730
https://github.com/HSU-ANT/ACME.jl/blob/main/src/ACME.jl#L733

So a couple of questions:

Q: Given the use of BigInt, are these tolerance filtering stages strictly neccesarily and/or is it because you'll run into problems when using the solution matrices and converting to floating point arithmetic later?

Q: (possible rephrasing of the above) - can we technically solve for the solution set purely symbolically, i.e. without component values for the matrices $D_v$, $D_i$ etc in:

(I'm already doing so for the linear/unique solution case, so I might expect this also to be the case for the nonlinear/under-determined case).

FInally in [2], there is the quote, which is where I'm getting the idea that tolerances might not be needed (due to exact arithmetic):

The proposed method has been implemented as part of the ACME framework. One of the main challenges encountered is the condition Fq,m,n = 0 which, when using floating point arithmetic, will usually only be fulfilled approximately. The obvious approach then is to treat entries with very small absolute value as zero. However, determining an appropriate threshold proves to be anything but trivial, as the scale of the entries in both q and z can differ by orders of magnitude from each other. So instead, we employ exact arithmetic, using rationals of arbitrary precision integers. This is possible as the whole model derivation process only needs a bounded number of additions, subtractions, multiplications, and divisions, so that the numerators and denominators may grow very large, but are still bounded.

I'm still getting my head around it all so apologies if this is phrased badly or doesn't make sense!

[1] "A generalized method for the derivation of non-linear state-space models from circuit schematics"
[2] "Automatic Decomposition of Non-linear Equation Systems in Audio Effect Circuit Simulation"

[martinholters]

I think you're pretty much on the right track: looking at this again, I agree that for Rational{BigInt}, one could just use tol=0 (and accordingly simplify those checks a bit). Note that float(Rational{BigInt}) == BigFloat and eps of that is about 1.73e-77, so it is quite unlike to stumble upon a value between 0 and tol in practice anyway.

And yes, this can absolutely be done symbolically. Instead of re-implementing gensolve, you could also implement Gaussian elimination to produce the reduced row echelon from and proceed from there. The advantage of gensolve is that it tries to exploit and preserve sparsity which probably also helps to keep the term size in check when working symbolically. Actually, it would be nice if gensolve just worked with Symbolics.jl, I might look into making that work.

[martinholters]

Actually, it would be nice if gensolve just worked with Symbolics.jl, I might look into making that work.

Hm, Symbolics and SparseArrays don't get along well, which is somewhat a showstopper here:

julia> using SparseArrays, Symbolics

julia> @variables x;

julia> sparse([x])' * sparse([1 0])
ERROR: TypeError: non-boolean (Num) used in boolean context

Interestingly, it works if the second argument is a sparse vector, not a sparse matrix, but I still don't want to rewrite gensolve to work around such issues.

A dedicated symbolic version that just does not use/work for SparseArrays shouldn't be too hard to do, but presently doesn't fit into ACME.

[martinholters]

So this is a version that seems to do ok for Symbolics:

function gensolve(a, b, x, h)
    m = size(a, 1)
    t = sortperm([count(!iszero, aᵢᵀ) for aᵢᵀ in eachrow(a)]) # row indexes in ascending order of nnz
    for i in 1:m
        aᵢᵀ = a[t[i], :]' # aᵢᵀ is a row of the a matrix
        s = aᵢᵀ * h
        # find the non-zero values in s
        jnz = findall(!iszero, s')
        if isempty(jnz)
            # aᵢᵀ is linearly dependent on the previous rows -> ignore it
            # note: if b[t[i],:]' ≠ aᵢᵀ*x, the equation system is inconsistent
            continue
        end
        # of those, pick the one where the column in h has fewest non-zero entries
        j = argmin(j -> count(!iszero, h[:,j]), jnz)
        q = h[:,j] / s[j]
        x = x + q * (b[t[i],:]' - aᵢᵀ*x)
        if size(h, 2) > 1
            h = h[:,[1:j-1;j+1:end]] - q * s[[1:j-1;j+1:end]]'
        else
            h = similar(h, (size(h, 1), 0))
        end
    end
    return x, h
end

gensolve(a, b) =
    gensolve(a, b, zeros(promote_type(eltype(a), eltype(b)), size(a, 2), size(b, 2)), Matrix{eltype(a)}(I, size(a,2), size(a,2)))

Not tested much, but for the one case I tried, it did work.

[hemmer]

Thanks so much for the response! Yes, the logic on line 730 never activates in practice, but with thresh=0.1 (default) there is filtering on line 733. I'll play with thresh=0 and assuming it all works fine will proceed from there re-implementing gensolve. I've not used Julia before (loving it though) so getting my head round that too.

My linear algebra is a little rusty so will need to revise that, but ultimately it'd be nice to switch to sympy built-ins etc once I've got a handle on gensolve. If you don't mind I'll keep the thread open to share any findings?

Appendix
My simple nonlinear problem to test my implementation with

using ACME
using Plots
fs = 44100

function saw(t; freq=10, ppv=10)
    x = @. ppv/2*(2*mod(t*freq, 1)-1)
    return reshape(x, 1, length(x))    
end;

function diode_filter(::Type{Circuit}; R=1000,C=9.400000000000001e-08)
    @circuit begin
    # elements
    j0 = voltagesource()
    r = resistor(R)
    c = capacitor(C)
    d1 = diode()
    jout = voltageprobe()

    # node connections
    #= node 00 =# j0[-] == c[2] == d1[-] == jout[-] == gnd
    #= node 01 =# j0[+] == r[1]
    #= node 02 =# r[2] == c[1] == d1[+] == jout[+]

    end
end

circuit = diode_filter(Circuit)
model = DiscreteModel(circuit, 1//fs)

ts = [n/fs for n in 0:10000]
x = saw(ts, freq=25)
y = run!(model, x)
print("Running diode_filter complete!\n")

plot(ts, x[1, :], label="V_in", size=(750,350))
plot!(ts, y[1, :], label="V_out", size=(750,350))

[martinholters]

Good point about

ACME.jl/src/ACME.jl

Line 733 in 57a512c

jat = jnz[nz_abs_vals .≥ thresh*max_abs_val] # cols above threshold

. Here thresh tries to strike a balance between sparsity promotion and numerical stability. Obviously, with arbitrary precision or symbolic arithmetic, numerical stability is not an issue, so yes, setting thresh=0 should be sound.

[martinholters]

I've written a brief summary of the algorithm together with a not-so-brief example of its application. That might be a better starting point for implementing it in another environment than the Julia implementation here.

[hemmer]

That's great, and super useful to have a reference solution.

If I enforce the same choice of J_n, then I recover your results (by default I get a different presumably equally valid solution, but that's also expected). And the method works fine for linear cases which is a good sanity check!

---

Would you be interested to continue this worked example for the most simple / naive solution method? I appreciate that in practice you want to apply dimensionality reduction (reduce_pdims!?) and other tricks (HomotopySolver, seperating and tracking multiple nonlinearities etc), but it makes the ACME.jl code a little hard to follow.

Even just a little more exposition/pseudocode of the high level maths might be useful for newcomers like myself, e.g. this is a sketch of what the Newton loop looks like, then you update q etc etc. I think it's because the argument q to the nonlinear function f(q) is not actually what you're iterating in the Newton steps I think (do you somehow iterate z instead?). There's also a reference to p in [1] but I'm not 100% sure on the relationship between z, q, p. Could even just be a single diode for simplicity.

No worries if not, I'm sure I'll work it out eventually!

[1] A GENERALIZED METHOD FOR THE DERIVATION OF NON-LINEAR STATE-SPACE MODELS FROM CIRCUIT SCHEMATICS