[SciPy-User] 1-norm estimation via Hager's algorithm in SCIPY

[Pauli Virtanen]

11.04.2015, 04:07, Wenlei Xie kirjoitti:
> I would like to estimate the 1-norm of matrix A*B, where A is n*k and B is
> k*n.  n is around 1 million while k is around 100. So I cannot materialize
> it directly.
> 
> I am wondering if there is anything similar to condest in Matlab to
> estimate the matrix 1-norm via Hager's algorithm? Would
> scipy.sparse.linalg.onenormest be a good choice? (It seems to be working
> for Sparse matrices?)

You can realize the matrix C=A*B as an abstract linear operator like this:

from scipy.sparse.linalg import aslinearoperator, onenormest
C = aslinearoperator(A) * aslinearoperator(B)
print(onenormest(C))

The resulting C however has a transpose defined only in Scipy versions
>= 0.15, so it won't be a valid input for onenormest in earlier versions.

Unfortunately, it seems that one step of the current implementation of
onenormest has quite large memory usage --- although it still scales
linearly with matrix dimension, there's a too big constant in front, so
it won't work OK for 100e6x100e6.