2
$\begingroup$

Let's say we want to solve for the eigen-values of a symmetric matrix of size $n$ x $n$.

In the Phase 1 of the computation, the matrix is reduced to a tridigonal form using Householder/Arnoldi's reduction.The theory says that this is $O(n^3)$.

In the Phase 2 of the computation, this tridigonal form is transformed to a diagonal form by using the Normalized Simultaneous Iteration (NSI). The two main steps in this iteraion are matrix multiplication of tridigonal matrices and the qr factorization of a tridigonal matrix. Each of these (if implemented specially for tridigonal systems) will be $O(n)$. Since these are to be repeated $n$ times, the total would be $O(n^2)$.

I want to see experimentally that the phase 1 actually reduces the time for computation. The easist way to do this is in MATLAB which uses the DGEMM routine for matrix multiplication and Householder reflector for qr factorization.

My question is : Does the qr algorithm and the DGEMM used in MATLAB take into account if the input matrix is tridigonal and optimize accordingly ? The documentation says that it differentiates between sparse and dense matrices.

If I have to write a code separately for matrix multiplication of tridigonal systems and a qr using householder reflector, will it be better than the ones provided in MATLAB ?

Is there a way to do this using something other than MATLAB , like a C or C++ library ?

$\endgroup$
  • $\begingroup$ Are you looking for routines that can efficiently compute eigenvalues of tridiagonal symmetric matrices? If so, have a look at dsteqr or dsterf (both are LAPACK routines). Both have an asymptotic ${\cal O}(n^2)$ cost. More information here. $\endgroup$ – GoHokies Nov 18 '15 at 18:58
  • $\begingroup$ Also, MATLAB has the built-in isbanded function that can detect the bandwidth of a given (sparse) matrix. It's not inconceivable to assume that they test for tridiagonal matrices somewhere in their qr implementation... $\endgroup$ – GoHokies Nov 18 '15 at 19:02
  • 1
    $\begingroup$ Have you tried benchmarking Matlab's eig, qr and *? Are they $O(n^3)$ or $O(n^2)$, with a tridiagonal matrix as input? It shouldn't take long to code. $\endgroup$ – Federico Poloni Nov 19 '15 at 5:58
4
$\begingroup$

Here's a short Matlab benchmark. I've used a sparse tridiagonal matrix.

function benchmark_qr

maxit = 25;
time = zeros(maxit,2);

for iter = 1:maxit

    mdim = 500 * iter;

    A = gallery('tridiag',mdim,-2,3,-2);

    t_start_qr = tic;
    qr(A);
    t_elapsed_qr = toc(t_start_qr);

    t_start_eig = tic;
    eig(A);
    t_elapsed_eig = toc(t_start_eig);

    time(iter,1) = mdim;
    time(iter,2) = t_elapsed_qr;
    time(iter,3) = t_elapsed_eig;

end

figure
loglog(time(:,1),time(:,2),'rs--')
hold on
loglog(time(:,1),time(:,3),'kd--')
loglog(time(:,1),(time(:,1)/100).^2,'--')
loglog(time(:,1),(time(:,1)/100).^3,'-.')
legend('qr','eig','N^2','N^3','Location','northwest');

end

The results suggest that both eig and qr are ${\cal O}(n^2)$ for sufficiently large $n$:

enter image description here

$\endgroup$
  • $\begingroup$ Why the loglog? The flop count is $f = {\cal O}(n^2)$, that is $f \le C n^2$ for sufficiently large $n$ (with $C$ a constant independent of $n$). Taking logarithms, we get $\log f \sim 2 \log n + \log C$ - a linear relationship with slope $2$ - something much easier to visualize and understand than the original (quadratic) equation. $\endgroup$ – GoHokies Nov 19 '15 at 9:30
  • $\begingroup$ Thanks. I got it. Now I'm trying the same thing for general matrices. I want to know if it is $O(n^3)$ or $O(n^2)$. $\endgroup$ – Srinivas K Nov 19 '15 at 10:12
  • $\begingroup$ qr is ${\cal O}(n^3)$ for general matrices. Define $A$ as A = rand(n) in my example to convince yourself that is the case. $\endgroup$ – GoHokies Nov 19 '15 at 11:27
  • $\begingroup$ Thanks a lot. As you said it looks like qr and even matrix multiplication is $O(n^3)$ for eneral matrices in MATLAB. $\endgroup$ – Srinivas K Nov 19 '15 at 14:01
  • $\begingroup$ @SrinivasK That should not be a surprise. $\endgroup$ – Jeff Nov 22 '15 at 1:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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