Given a large square matrix A and a sufficiently regular function f so that f(A) is well defined, we are interested in the approximation of the leading singular values and corresponding left and right singular vectors of f(A), and in particular in the approximation of ‖f(A)‖, where ‖⋅‖ is the matrix norm induced by the Euclidean vector norm. Since neither f(A) nor f(A)v can be computed exactly, we introduce a new inexact Golub–Kahan–Lanczos bidiagonalization procedure, where the inexactness is related to the inaccuracy of the operations f(A)v, f(A)⁎v. Particular outer and inner stopping criteria are devised so as to cope with the lack of a true residual. Numerical experiments with the new algorithm on typical application problems are reported.