lucaferranti
commented
2 years ago thank you for your interest in the project and sorry for the long delay in answering, very busy week :D
I think your idea sounds very interesting and would make sense to include something like that in the package (maybe the names of the API can be negotiated :) ). I'll try here to summarize your proposal to make sure I understand correctly and ask a few clarifying questions.
If I understand correctly you are proposing that given a matrix M, find a non-interval matrix P and an interval diagonal matrix D so that M is contained in P*D*P^{-1}. I think this would be interesting to explore and maybe it might have some practical utilities e.g. to bound functions of matrices.
-
Do I understand correctly, that the advantage of your proposed approach is that you are relaxing the problem as you are not requiring
Pto bound the eigenvectors, computing this factorization would be more efficient thanverify_eigenor the approach proposed in the paper in #116 ? I have a couple of scenarios in mind where this could be useful -
I could not understand from your description whether the matrix
Mwould be a matrix of floats (single points) or intervals. Which one did you have in mind? These would have slightly different problem formulations:- find
P, Dso thatM ∈ PDP⁻¹(starting from floating point) - find
P, Dso thatM ⊆ PDP⁻¹(starting from interval) - The reason I'm asking is because while the latter is more generals and algorithms working for that work also for the former, if you can assume that the starting point is a matrix of floats (or very small intervals), then you can generally find more performant algorithms, so I think it could be good to separate the two cases and start benchmarking only with the former could be an easier start.
- find
-
In your proposal you mention to compute
P⁻¹, as you probably know matrix inversion cannot be computed exactly with floating point, if you want to work with intervals and want rigorous computations generally there are two options:- Compute
P⁻¹rigorously, meaning in the factorizationPwould be a matrix of floats,Dan interval matrix andP⁻¹an interval matrix guaranteed to contain the true inverse ofP - Slightly modify the problem: Find
P, D, Rso thatM ∈ PDR, whereDis interval andP, Rnot. Then probably an approximate inverse ofPwould be good enough.
- Compute
-
Did you have in mind an algorithm on how to determine
PandD? Can you link some papers I should link to better understand this issue?
If you have further questions / comments, do not hesitate to ping me.
Feature description
I'm new to interval computing. I'm a PhD student at Imperial College London. I have a suggestion, and I'm wondering how you handle this.
Certain operations on matrices are uncomputable. For instance, eigendecomposition of even "nice" matrices like symmetric matrices is not computable. What "uncomputable" means in practice is that "forwards numerical stability" is not attainable. The notion of "backwards numerical stability" may still be attainable and can suffice for certain use cases. Unfortunately, interval arithmetic in the naive sense cannot "handle" backwards stability. But this limitation can be by-passed.
My suggestion surrounds eigendecomposition. The API I propose should have three functions:
eigendecomp : Mat(Complex) -> (Mat(Complex), DiagMat(Complex)),inv_eigendecomp : (Mat(Complex), DiagMat(Complex)) -> Mat(Complex)lift : (Complex -> Complex) -> ((Mat(Complex), DiagMat(Complex)) -> (Mat(Complex), DiagMat(Complex))).eigendecomp(M)should produce a pair of matrices(P,D)such that the matrix P is exact (and not interval valued!!) and D is interval-valued. P and D should be chosen so thatP D P^-1containsMas snuggly as possible. The actual eigenvalues and eigenvectors of M don't matter because the eigenvectors of M in particular are not computable.I also suggest that the purpose of eigendecomposition is to lift functions of type
Complex -> Complexto complex matrices. I suspect that the eigenvalues and eigenvectors are perhaps not entirely relevant, in practice.