# Rigorous spectral opnorm

**URL:** <https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739>\
**Category:** General Usage\
**Tags:** question, intervals, linearalgebra\
**Created:** [June 14, 2022, 11:18am UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739 "2022-06-14T11:18:16Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 14, 2022, 11:18am UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/1 "2022-06-14T11:18:16Z")

</div>

I am wondering if there exists a way to rigorously compute spectral norms (an upper bound) of matrices (so the usual 2-norm for matrices) using the IntervalArithmetics package ([GitHub - JuliaIntervals/IntervalArithmetic.jl: Rigorous floating-point calculations using interval arithmetic in Julia](https://github.com/JuliaIntervals/IntervalArithmetic.jl)). The matrix I am considering is a matrix of intervals.  
IntervalMatrices.jl allows the computation of the L^1 and L^\infty norms (as Steward pointed out) but not the spectal one.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 14, 2022, 1:15pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/2 "2022-06-14T13:15:38Z")

</div>

> [@mathmath](#):
>
> The LinearAlgebra package does not seem to use Interval computations.

The LinearAlgebra package computes things like `opnorm(A)` up to floating-point error. You only need interval arithmetic if you need bounds on the roundoff error of the computation.

To get bounds on the roundoff errors, you need interval-arithmetic algorithms to bound the singular values of A (or the eigenvalues of A^\* A). Such algorithms exist — see e.g. [Hladik (2009)](https://hal.inria.fr/inria-00370603/document), [Hladik (2013)](https://doi.org/10.1016/j.amc.2012.11.075), [Nerantzis & Adjiman (2015)](https://doi.org/10.1080/10556788.2016.1184663), [Hartman et al. (2019)](https://arxiv.org/pdf/1912.05275.pdf), and probably others — and the Hladik (2013) algorithm is [implemented in Julia](https://juliaintervals.github.io/IntervalLinearAlgebra.jl/dev/api/eigenvalues/). Moreover, it has a function `bound_perron_frobenius_eigenvalue` to bound just the maximal eigenvalue of an SPD matrix like `A'A`, which is precisely what you need for the operator norm. I’m not sure how tight of a bound it returns, however — the usual caveats about interval arithmetic potentially returning very loose bounds applies.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 14, 2022, 1:20pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/3 "2022-06-14T13:20:18Z")

</div>

It seems that the [IntervalMatrices.jl package](https://juliareach.github.io/IntervalMatrices.jl/latest/lib/methods/#Norms) already provides an `opnorm` function to bound the matrix norm of an interval matrix, but apparently [it currently only supports](https://github.com/JuliaReach/IntervalMatrices.jl/blob/6d130cf83e946375ed2ac8748b0d4c15ad297f16/src/operations/norm.jl#L21-L29) the induced (operator) L^1 or L^\infty norms.

---

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 14, 2022, 1:40pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/4 "2022-06-14T13:40:55Z")

</div>

Thank you for the reply, indeed I need a rigorous bound on the roundoff error.  
I have two questions :

- the algorithms you gave enables to verify an eigenvalue of A\*A, but it does not verify that it is actually the biggest one ?
- for the Perron Frobenius eigenvalue, doesn’t it work only for elementwise non-negative matrices ?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 14, 2022, 2:16pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/5 "2022-06-14T14:16:53Z")

</div>

> [@mathmath](#):
>
> - the algorithms you gave enables to verify an eigenvalue of A\*A, but it does not verify that it is actually the biggest one ?

`eigenbox` apparently encloses _all_ the eigenvalues of a matrix, including the biggest one.

> [@mathmath](#):
>
> for the Perron Frobenius eigenvalue, doesn’t it work only for elementwise non-negative matrices ?

Yes, sorry, this won’t work for general SPD A^\*A.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 14, 2022, 2:18pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/6 "2022-06-14T14:18:47Z")

</div>

> [@mathmath](#):
>
> Thank you for the reply, indeed I need a rigorous bound on the roundoff error.

Do you need a lower bound on the operator norm, an upper bound, or both? A lower bound is _much_ easier to compute.

---

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 14, 2022, 2:35pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/7 "2022-06-14T14:35:32Z")

</div>

The problem is that I am treating huge matrices (10 000 by 10 000) so computing all the eigenvalues will be too much in term of computations. That’s why I am trying to find out a method that only focus on the maximum eigenvalue (or spectral radius) to avoid unecessary computations.

---

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 14, 2022, 2:38pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/8 "2022-06-14T14:38:22Z")

</div>

I need on upper bound for the operator norm. But as my matrix is invertible I can also compute the minimal eigenvalue of the inverse. Then, a lower bound for such an eigenvalue will give me an upper bound for the original problem. I don’t know if this helps.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 14, 2022, 2:52pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/9 "2022-06-14T14:52:16Z")

</div>

Do you really need super rigorous bounds on the opnorm of a 10k by 10k matrix? `opnorm` is already super accurate. Presumably if your system is of size 10k you have much bigger problems than that (eg discretization issues)

---

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 14, 2022, 2:58pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/10 "2022-06-14T14:58:22Z")

</div>

The thing is that I need the bound to be rigorous. So I do not necessarily need it to be super sharp (even if that would be better 🙂 ) but simply rigorous.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 14, 2022, 3:11pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/11 "2022-06-14T15:11:37Z")

</div>

One thing you can do to verify a bound is making sure `A'A-lambda I` is negative. Don’t know if that helps, but at least it reduces the problem to checking whether an interval symmetric matrix is positive, which might be a bit faster? But honestly, life is much easier when you embrace floating point error 🙂

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [June 14, 2022, 8:08pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/12 "2022-06-14T20:08:35Z")

</div>

You could try a hybrid approach: just compute the SVD in the normal way and then write out a perturbation bound for the largest singular value. The perturbation theory is going to involve quantities (e.g. residuals) that will involve further computation involving A, the singular value, and its singular vectors. You might then worry about rounding errors in that computation undermining the rigor of your bounds. However you could apply interval arithmetic to computing the residuals and quantities in the perturbation bound, which would be quite a bit nicer than applying it to the SVD. I haven’t worked seriously with interval arithmetic before, but I did have the impression that the intervals could grow exponentially in computing matrix decompositions. Since the largest singular value is perfectly well conditioned, the perturbation theory reflects that, and the residuals etc. won’t be too bad in interval arithmetic, you would probably get a much better bound relative to trying to apply interval arithmetic directly to the SVD. I don’t mind sketching it out if you want to give it a try.

I hope that makes sense. I think it should work well, but I haven’t tried anything like that before and I have more experience with perturbation bounds than interval arithmetic.

Edit: This doesn’t work if you are starting with an interval matrix. But I had the impression from your wording that you might be starting with a matrix of numbers instead of intervals.

---

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 14, 2022, 8:31pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/13 "2022-06-14T20:31:26Z")

</div>

In fact I am starting with an interval matrix, my description was not clear…

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [June 14, 2022, 8:57pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/14 "2022-06-14T20:57:29Z")

</div>

Hmmm… I guess I was an optimist and headed for the easier problem. A follow-up question would be: Is this an interval matrix that represents small intervals added to a larger matrix? (i.e. a fixed matrix plus smaller intervals representing a smaller perturbation) Or is it really intervals that are larger?

---

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 14, 2022, 11:07pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/15 "2022-06-14T23:07:19Z")

</div>

Exactly, the intervals are pretty small, in practice they are of the order of the machine precision. I am building my matrix in Floats then turning the entries in Intervals before computing the norm, so it is as good as it can get.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [June 15, 2022, 11:24am UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/16 "2022-06-15T11:24:35Z")

</div>

I think the details are important here. So, just to make sure, here is my understanding: You have an interval matrix {\bf A} = [\underline{A}, \overline{A}] and you want an interval bound on \|A\|\_2 that holds for all A\in {\bf A}? Further, \overline{A} - \underline{A} is of the order of the machine precision. So any computations you do in computing the final interval for \|A\|\_2 are going to introduce errors that are comparable to the width of the interval you are trying to compute. So those must be taken into account as well? I’m using notation that I think is fairly standard and is summarized in the 2009 paper by Hladik that @stevengj cited. And the problem is large enough that you want to look at the largest singular value alone and don’t consider computing a full SVD to be a reasonable step in the process?

If you ignore the additional numerical errors and just want a decent upper bound on \|A\|\_2, I think `sqrt(eigenbox(Symmetric(AI'*AI)))` for an interval matrix `AI` will give that. But they are explicit in saying that “The enclosure is not rigorous, meaning that the real eigenvalue problems solved internally utilize normal floating point computations.” It wasn’t that slow on my very old system. It seemed a little memory hungry, so I couldn’t try 10,000\times 10,000, but it took about 3 minutes for a 7,000\times 7,000.

If you can deal with a couple of SVDs, I think the perturbation theory approach could handle both the intervals and the numerical errors. Maybe you could apply the same thing to a singular values computed using a power method, but it is an additional complication in something that is already a bit tricky.

---

<div class="post-metadata">

**Author:** ![mathmath](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mathmath/32/37157_2.png) [@mathmath](https://discourse.julialang.org/u/mathmath)\
**Post date:** [June 15, 2022, 12:00pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/17 "2022-06-15T12:00:11Z")

</div>

Your description is exactly what my problem is about.

> So any computations you do in computing the final interval for ||A||\_2  
> are going to introduce errors that are comparable to the width of the interval you are trying to compute. So those must be taken into account as well?

Yes indeed, the standard functions (sum, multiplication,…) are already implemented to be rigorous in IntervalMatrices or IntervalLinearAlgebra packages but at some point the width of the intervals will come into place.

> And the problem is large enough that you want to look at the largest singular value alone and don’t consider computing a full SVD to be a reasonable step in the process?

Yes, computing a full rigorous SVD will be to much in term of computations (also is there a package that enables the computation of a rigorous SVD ? the IntervalLinearAlgebra package does not seem to do so). Computing a SVD on a block of the matrix and trying to bound the rest could be an option, but that is still a complicated problem.

> If you can deal with a couple of SVDs, I think the perturbation theory approach could handle both the intervals and the numerical errors.

What do you mean by dealing with a couple of SVDs ?

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [June 15, 2022, 12:28pm UTC](https://discourse.julialang.org/t/rigorous-spectral-opnorm/82739/18 "2022-06-15T12:28:55Z")

</div>

If you can afford to compute a couple of ordinary singular value decompositions of (non-interval) matrices A\_c = (\overline{A} + \underline{A})/2 and A\_\Delta = (\overline{A} - \underline{A})/2, I think you can apply perturbation theory in a way that gives bounds that handle both the numerical errors and the original intervals. A couple of 10,000\times 10,000 SVDs shouldn’t be a deal breaker on modern computers, unless you are doing this lots of times. (Which maybe you are, since you are so worried about performance…)

You could apply similar perturbation theory to just the maximum singular values for A\_c and A\_\Delta computed using some other method (e.g. the power method or Lanczos.) But you are then going to have to worry about whether you converged to the correct singular values. It’s hard to compete with a full matrix decomposition if you want to be absolutely sure. If the second singular value is close to the first and has a slightly larger interval, maybe you really wanted the second singular value anyway. I think you are almost inevitably going to sacrifice some certainty if you aren’t looking at the smaller singular values. If the second singular values are known to be significantly smaller than the first, you might get something that is a strict bound with high probability, barring some terrible choice of starting vector for your iteration.