MATLAB has a couple of "exact" functions for this, cond
and rcond
, with the latter returning a reciprocal of the condition number. Matlab approximate function condest
is more fully described below.
Often estimates of the condition number are generated as by-products of the solution of a linear system for the matrix, so you might be able to piggyback the condition number estimates on other work you need to do anyway. See here for a brief description of how estimates are computed. Also Sandia Labs AztecOO documentation remarks (see Sec. 3.1) that optional condition number estimates are available from iterative solvers (using the generated tridiagonal Lanczos matrix with Conjugate Gradients or the generated Hessenburg matrix with Restarted GMRES).
Since your matrices are "very large" and "only available as functions", the logical approach would be a method that piggybacks on a conjugate gradient solver or variant.
A recent arXiv.org paper Non-stationary extremal eigenvalue approximations in iterative
solutions of linear systems and estimators for relative error proposes such an approach and has a few citations to the earlier literature.
Now that I look, this forum has a number of closely related previous Questions (not all with Answers, but check Comments):
Estimate extreme eigenvalues with CG
Estimation of condition numbers for very large matrices
Fastest algorithm to compute the condition number of a large matrix in Matlab/Octave
Since availability of MATLAB code was part of the Question, here's some information about condest
, a built-in function that estimates the 1-norm condition number $\|A\|_1 \|A^{-1}\|_1$. The idea is from Hager(1984), with a 2010 write-up and extensions here, to explicitly compute $\|A\|_1$ (find maximum 1-norm of a column) and estimate $\|A^{-1}\|_1$ by a gradient method. See also John Burkardt's CONDITION, a MATLAB library (other languages available) "for computing or estimating the condition number of a matrix."
Since your matrix is apparently Hermitian and positive definite, perhaps the 2-norm condition number is of greater interest. The problem then amounts to estimating the ratio of largest to smallest (absolute) eigenvalues. The challenge is somewhat parallel to the 1-norm case in that generally a good estimate for the largest eigenvalue can be easily obtained, but estimating the smallest eigenvalue proves more difficult.
Although aiming at non-SPD (and even non-square) cases, this recent arXiv.org paper, Reliable Iterative Condition-Number Estimation, gives a good overview of the smallest eigenvalue estimation problem and a promising line of attack by a Krylov-subspace method (LSQR) that amounts to Conjugate Gradients in the SPD case.