15

The Mathematica documentation says it is possible to estimate the matrix condition number in norms 1, 2, and ∞. But the 2-norm raises a message.

This is an extract from reference documentation "tutorial/LinearAlgebraMatrixComputations" where I changed the expression to compute the 2-Norm.

error

UPDATE

As requested, if you want to try by yourself and check the error on your system, simply type:

mat = {{1., 2.}, {3., 4.}};
LinearAlgebra`MatrixConditionNumber[mat, Norm -> 2]
Community
  • 1
unlikely
  • 7,103
  • 20
  • 52
  • 2
    You can compute it yourself with Norm[mat,2]Norm[Inverse@mat,2], but this is unsatisfying... – jtbandes Jul 17 '14 at 06:27
  • @jtbandes yes, possible, but unsatisfying... – unlikely Jul 17 '14 at 06:36
  • You can estimate the L_infinity condition number of a square matrix via LUDecomposition (the third part of the result is that estimate). For the 1-norm take the LUDecomposition of the transpose. – Daniel Lichtblau Jul 17 '14 at 21:47

3 Answers3

14

The compatibility information at Compatibility/tutorial/LinearAlgebra/MatrixManipulation says

These functions were available in previous versions of Mathematica and are now available on the web at library.wolfram.com/infocenter/MathSource/6770:

LinearEquationsToMatrices
InverseMatrixNorm
ConditionNumber

You can download the original package there. It's too long to provide an excerpt here, but you can load it and use it in your code as-is.

(There seems to be a vestigial version of LinearAlgebra`MatrixConditionNumber which, as you noticed, only supports norms 1 and ∞.)

On the other hand, if you are okay with the computation involved in producing an exact answer, the documentation for SingularValueList says

The 2-norm of a matrix is equal to the largest singular value … The 2-norm of the inverse is equal to the reciprocal of the smallest singular value … [Thus,] The condition number of the matrix is equal to the ratio of largest to smallest singular values.

So you can use:

First@#/Last@#& @ SingularValueList[mat]


Performance

The old implementation for 2-norms is considerably faster for large random matrices (it is worth noting that both implementations seem to take advantage of multiple cores):

performance

While the relative error stays low (this may depend on the precision of your input):

error

jtbandes
  • 1,422
  • 1
  • 11
  • 20
  • I'm not against an exact answer, but I'm against the time required for an exact answer if the matrix is very large. The documentation suggest the MatrixConditionNumber function returns a fast, approximate, answer, for large, poterntially banded or sparse, matrices, and this is what I'm searching for. I'll give a try to the old package... – unlikely Jul 17 '14 at 06:53
  • The last edit is detailed and interesting. A question that arise is why this old function is now deprecated (and also why the documentation is out of sync): maybe there are some issue using this estimate? For example I noticed I often get large errors in the estimate with small matrices. And for small matrices the exact computation is often faster. – unlikely Jul 17 '14 at 12:32
  • The best idea is probably to extract from the old package the relevant definitions for the estimate of the spectral condition number, and write a helper function that dispatch the request to the built-in Norm (for 1 and Infinity norm), to the old package function for the 2-Norm and large matrices, to the exact computation function for the 2-Norm and small matrices... Strange that Wolfram forgot to do this... – unlikely Jul 17 '14 at 12:33
  • That sounds reasonable; the threshold would of course be up to you, based on the computational speed you need. – jtbandes Jul 17 '14 at 16:41
  • For complete accuracy, you need to set a particular option. Also, since you don't need all the singular values, you can just generate the smallest and the largest separately: cond2[mat_?MatrixQ, opts___] := SingularValueList[mat, 1, opts]/SingularValueList[mat, -1, Tolerance -> 0, opts] – J. M.'s missing motivation Jul 22 '15 at 14:40
8

Other solutions are fine, but they use old MatrixConditionNumber as a magic box. However, it has a simple idea. The 2-norm condition number of the matrix $M$ is a ratio $\sigma_{\rm max}/\sigma_{\rm min}$ between the maximum and the minimum singular values. The maximum singular value $\sigma_{\rm max}$ can be estimated by a simple power iteration:

\begin{align} u_0&={\rm random},\\ u_{i+1}&=M M^\dagger u_i/|u_i|,\quad i=1,\ldots,n-1\\ \sigma_{\rm max}^2&\approx |u_n|. \end{align}

The minimum singular value $\sigma_{\rm min}$ can be estimated by the same power iteration with $M^{-1}$ instead of $M$, which means to solve linear system on each step. So

condNum2[m_, k_: 10] := Module[{s, u1, u2},
  s = Internal`DeactivateMessages@LinearSolve@m;
  {u1, u2} = RandomReal[NormalDistribution[], {2, Length@m}];
  Do[
   u1 = m.Conjugate[Conjugate@Normalize@u1.m];
   u2 = s[s@Normalize@u2, "C"];
   , {k}];
  Sqrt[Norm@u1 Norm@u2]
  ]

SeedRandom[0];
n = 1000;
m = SparseArray[RandomInteger[{1, n}, {10 n, 2}] -> RandomComplex[1 + I, 10 n]];

condNum2[m] // AbsoluteTiming
(* {0.205369, 3599.44} *)

MatrixConditionNumber[m, 2] // AbsoluteTiming
(* {0.220560, 3599.45} *)

Here s[...,"C"] solves the conjugated transposed linear system. Another possibility is to find the maximum and the minimum singular values by the Arnoldi algorithm 

SingularValueList[m, 1, Method -> "Arnoldi"]/
   SingularValueList[m, -1, Method -> {"Arnoldi", "Shift" -> 0}] // First // AbsoluteTiming
(* {0.464198, 3599.46} *)
ybeltukov
  • 43,673
  • 5
  • 108
  • 212
  • The complexity is still that of solving linear systems so it involves at minimum something like LUDecomposition for dense matrices and perhaps an iterative method for sparse. – Daniel Lichtblau Jan 31 '15 at 21:18
  • 1
    I wouldn't recommend the power iteration in general, at least in the form given here. Forming the Gram matrix, even implicitly, squares the condition number, which defeats the purpose of why you're looking at the condition number of the original in the first place. A concrete example would be the Läuchli matrix. – J. M.'s missing motivation Jul 22 '15 at 14:28
6

The LinearAlgebra package has been deprecated since Mathematica 5, and is no longer bundled with Mathematica 9 or newer. You can still download a part of it (which contains the MatrixConditionNumber function) at the URL that jtbandes gave in his answer.

First, we need to load the package:

(* Be sure to install the above linked package in a "LinearAlgebra" subfolder for M9. *)
<< LinearAlgebra`MatrixManipulation`

Now if you call MatrixConditionNumber without a context in front it works and with the second argument as a single number instead of an option, it works:

MatrixConditionNumber[{{1., 2.}, {3., 4.}}, 2]

14.933

Note that you can do this without the LinearAlgebra package as follows:

Max[#]/Min[#]& @ SingularValueList[ Inverse[{{1., 2.}, {3., 4.}}] ]

14.933

Teake Nutma
  • 5,981
  • 1
  • 25
  • 49