I am going through your paper. I need the smallest singular value of the matrix to compute its 2-norm condition number. I cannot use SVD due to memory issues. I computed the largest singular value using the power method, which converged very fast. However, the inverse power method used to compute the smallest singular value is very slow and is not converging.

I have another matrix of the same dimension. From modeling perspective I know that this matrix is well conditioned. However, I need to quantify my presumptions. Inverse power method converged for this matrix but it took time.

I came to the conclusion that GMRES+ILU as a matrix inverter was not sufficient for the previous ill-conditioned matrix. Adding a good preconditioner may also take time.

Right now, I am going through implicitly restarted Lanczos bidiagonalisation algorithm. Since I am new to this field I do not have much knowledge about the state of the art techniques used for this type of computations.

How about using the Lanczos method for A^T A?

Please try the keyword "LOBPCG" (locally optimal block preconditioned conjugate gradient). The method was developed by A. Knyazev e.a. It is based on the Rayleigh quotient minimization with the use of  the preconditioned block CG method. Of course, the Rayleigh quotient (as well as the preconditioner) must be written for A^TA in your case.

Most remarkable, there must be some free software implementing this method.

did you mean ATA of the augmented form ?

What is "augmented form"? I mean the same thing as David E. Stewart :

reformulate your problem as the standard symmetric eigenproblem and use some standard software.

If you will provide some appropriate preconditioner for A^TA or AA^T, then LOBPCG type methods can be very efficient.

I don't know if you meant to find ATA directly. Directly finding ATA in my case is difficult due to the following reasons:

1) My matrix size is very large. Although A is sparse, ATA may become dense. This may lead to storage issues.

2) If the singular values are near to machine epsilon it may lead to inaccurate computation.

In order to circumvent these problems symmetric eigenvalue  computation is performed on the augmented form C(column vector)= (C1  C2)T where C1(row vector)= (0 A), C2(row vector)=(AT 0). Sorry for the above misleading formulation. C is a block matrix of size 2nx2n when size of A is nxn with two zero matrices through the left diagonal and the other two blocks are A and AT.

In the Lanczos method and other iterative methods, it is not necessary to form A^TA or AA^T explicitly. One only needs to perform matrix-vector multiplications such as z=A^T(Av). Or such an operation as z=A^Tw is also unavailable? Some preconditionings for A^TA can also be constructed without any explicit forming of A^TA.

Yes. LOPCG works. Thanks.

This approach is now implemented in https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html