3
$\begingroup$

enter image description hereI am trying to do SVD of a large block-hankel matrix for model order reduction (Low rank approximation). However, I quickly run into memory issues in forming the large Block-Hankel matrix and CPU issues due to running SVD itself.

In theory, it seems like we don't have to form the Hankel matrix itself. As long as we have the elements of either the first row/column of the matrix, we have complete information of all the entries of the Hankel matrix, and we should be able to perform the SVD without even forming this large matrix in memory ?

Is there any algorithm (in netlib or any of the other numerical libraries) or any other user-submitted MATLAB package that takes advantage of the special structure of the block-Hankel matrix ?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
4
$\begingroup$

There are matrix free algorithms (algorithms that use only matrix-vector multiplications rather than working directly with the entries of the matrix) that can compute approximate values of a few singular values/vectors of a large matrix. Since you want a low rank approximation to this matrix, you could use such an algorithm to find the $k$ largest singular values and associated singular vectors.

In order to make this work you will need to develop routines that can compute matrix-vector productions $Ax$ and $A^{T}y$ for your structured matrix $A$.

ARPACK is perhaps the most widely used package in this area. In MATLAB if you have $A$ as a sparse matrix you can use the svds() function which calls ARPACK routines. Unfortunately, MATLAB doesn't seem to provide an interface to these functions that allows you to just provide matrix-vector multiplication routines instead of a sparse matrix.

There's a survey of available software package for sparse eigenvalue and singular value problems that might be helpful to you. I assume that someone has put together MATLAB interfaces to some of these libraries that do have the ability to work with user supplied matrix-vector multiplication routines.

As I recall, PROPACK includes a MATLAB interface that allows for user supplied matrix-vector multiplication routines.

Improve this answer
$\endgroup$
2
  • $\begingroup$ This is indeed tremendously useful....Thank you ! However, since I am an electrochemist and not a linear algebra expert, I am a little lost when you say matrix-Vector product can be used to come up with SVD. Hence, I'd appreciate a little bit more low-level detail on how to actually code up these things. I simply call svd(A), & MATLAB produces [U,S,V] without worrying about details of the numeric implementation. I also understand that svds() in MATLAB is more helpful for sparse computations. Although I gather that your answers are helpful, i appreciate any further help. $\endgroup$
    – Dr Krishnakumar Gopalakrishnan
    Jul 11, 2015 at 10:41
  • 1
    $\begingroup$ PROPACK and similar packages use iterative algorithms in which some of the computational steps consist of multiplying the matrix A times a vector. (Also, it's necessary to multiply $A^{T}$ times a vector.) Since you don't have sufficient storage to hold the entire matrix, you can't simply use MATLAB's built-in operation to do Ax. Rather, you will have to write subroutine that that takes as input a vector x (of length 70,000) and outputs Ax (a vector of length 350,000) using the blocks of your block Hankel matrix. $\endgroup$
    – Brian Borchers
    Jul 11, 2015 at 17:07
0
$\begingroup$

I suggest to avoid the SVD, use Lanczos with Fast Block matrix-vector multiplication via FFTs. You don't need to save the matrix, it is all on the flight. The theory is explained here

A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions,Jianjun Gao1, Mauricio D. Sacchi and Xiaohong Chen, GEOPHYSICS,Volume 78, Issue 1. http://library.seg.org/doi/abs/10.1190/geo2012-0038.1

Cheers, Mauricio D Sacchi

Improve this answer
$\endgroup$
1
  • $\begingroup$ I suspect you are suggesting that for low-rank approximation, only the largest singular values are needed and can be found less expensively than the full decomposition. However your brief post leaves much to the Readers imagination. $\endgroup$
    – hardmath ♦
    Jul 25, 2015 at 3:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.