Computing the logarithmic spectral norm rapidly

I wish to compute the logarithmic spectral norm of a square matrix $A$, which is defined by

$$\mu_2[A]=\lim_{t\downarrow0}\frac{\|I+tA\|_2-1}{t}.$$

Hence, I simply code the above formula in Mathematica, but the deal is for large sizes, the computation will not stop and does not provide any results! I think, a kind of $\frac{0}{0}$ happens in the middle of the process and that is why this raw implementation does not work properly.

 ClearAll["Global*"]
logarithmicNorm[B_] :=
Limit[(Norm[Id + t*B, 2] - 1)/t, t -> 0., Assumptions -> t > 0]

n = 2;
Id = SparseArray[{{i_, i_} -> 1.}, {n, n}, 0];
A = N@SparseArray[{{i_, i_} -> 2.6, {i_, j_} /; Abs[i - j] == 1 ->
RandomReal[]}, {n, n}, 0];

logarithmicNorm[A]


The above sample example works fine when $n=2$, but if I choose $n\geq3$ it gets stuck! Now, the point is I want to compute this kind of norm for matrices of high dimensions, such as $n>200$.

I will be thankful if anyone gives me some tips or help.

• I'm pretty sure that Limit is a function that evaluates limits purely symbolically, which means it has to convert to a function first, then if it knows that functional form, it can take the limit. You will probably need to design your own limit function that numerically finds a large enough t such that the value of the function doesn't change very much (within some user-defined tolerance). Commented Jan 25, 2016 at 19:16

The formula for spectral norm you are using is meant to be the formal mathematical definition of the quantity. However this is restrictive for practical use as symbolic norm calculation on high dimensions are very cumbersome.

The formulation you might be looking for is the following. Here $\mu_{2}$ is the logarithmic two norm. $$\mu_{2}(A) := \lambda_{max}\left(\frac{A+A^{\mathrm T}}{2}\right)$$ where $\lambda_{max}(M) = \max_i\{ {\rm Re}(\lambda_i) \}$ is the spectral abscissa of a square matrix $M$.

So now you can do the following.

n = 1400;
SeedRandom[1234];
A = N@SparseArray[{{i_, i_} -> 2.6, {i_, j_} /; Abs[i - j] == 1 ->
RandomReal[]}, {n, n}, 0];
logarithmicNorm[A_] :=
Eigenvalues[(A + Transpose@A)/2, 1,
Method -> {"Arnoldi", "Criteria" -> "RealPart", "BasisSize" -> 400,
"MaxIterations" -> 10^5} ];
logarithmicNorm[A]


{4.35321}

You can increase the "BasisSize" and "MaxIterations" accordingly if you are working with large sparse matrices.

• Wow. It gives correct results, but may you please let us know where are "the limit" or "the identity matrix" (which are in the mathematical definition) are?
– Faz
Commented Jan 25, 2016 at 20:02
• How the limit is magically removed and a simplified formula takes the center-stage is out of scope of this answer. But it can be proven to be so using functional analysis apparatus on a Hilbert Space of operators. Here is the original paper by C. V. Pao. By now it is a standard result in the fields of dynamical systems and their stability. Commented Jan 25, 2016 at 20:43
• Great. Thank you so much.
– Faz
Commented Jan 26, 2016 at 11:33
• Can you elaborate why you chose these specific settings for Eigenvalues? Commented Nov 3, 2016 at 11:21
• As far as I can remember...just to highlight the options of Eigenvalues. Also for convergence I probably needed a higher "MaxIterations"` than the default one. Commented Nov 3, 2016 at 14:34