4
$\begingroup$

I am looking for a good numerical method to solve a linear system (of small/moderate size) with a matrix of the form $$ \begin{bmatrix} a_1 & b_1\\ & a_2 & \ddots \\ & & \ddots & b_{n-1}\\ b_n & & & a_n \end{bmatrix}. $$ What would you suggest?

Entries not displayed are zero. The $a_i$ and $b_i$ entries can be arbitrary, so I would not trust the stability of unpivoted Gaussian elimination.

Related question: solving system of linear equations with cyclic tridiagonal matrix.

$\endgroup$
  • 1
    $\begingroup$ Since you say you have a lot of these, would it be worth trying, for a given matrix $A$, to use unpivoted GE first, together with computing an estimate of its accuracy, and if it fails or it's too inaccurate, only then compute pivoted LU. This would take time $O(n+\alpha n^3)$, where $\alpha$ is the fraction of ill-behaved matrices, which, if small enough, makes it linear. $\endgroup$ – Kirill Oct 21 '16 at 18:27
2
$\begingroup$

The special structure of your matrix is easily exploited by a custom made $QR$ factorization based on Givens rotations. This leads to a method of $O(n)$ complexity. Below is a Octave/Matlab code, which solves $Mx=rhs$ for cyclic bidiagonal matrices $M$ of size $n \ge 2$ .

% Create a cyclic bidiagonal testmatrix
a=[2 -3 -1 5 3]';
b=[7 2 -4 -1 6]';
rhs=[-1 -3 2 6 2]';
n=5;

% Full matrix solution
M=diag(a)+diag(b(1:n-1),1);  M(n,1)=b(n);
x_full_matrix=M\rhs;
res_nrn_full_matrix=norm(M*x_full_matrix-rhs)

% QR factorization using Givens rotations (see Golub and Van Loan), without 
% forming the full Matrix M. M is transformed to an upper triangular 
% matrix R (which is also not formed as a full matrix)
% a and b are overwritten. d stores fill-in entries of R
c=zeros(n-1,1);
s=zeros(n-1,1);
d=zeros(n-2,1);

z=b(n);
for i=1:n-1,
  % Compute Givens rotation
  if z==0, c(i)=1; s(i)=0;
  else
    if abs(z)>abs(a(i)), tau=-a(i)/z; s(i)=1/sqrt(1+tau*tau); c(i)=s(i)*tau;
    else                 tau=-z/a(i); c(i)=1/sqrt(1+tau*tau); s(i)=c(i)*tau;
    end
  end
  % Apply Givens rotation
  a(i)=c(i)*a(i)-s(i)*z;
  a_n=a(n);
  if i==n-1, break; end
  z   =s(i)*b(i);
  b(i)=c(i)*b(i);
  a(n)=c(i)*a_n;
  d(i)=-s(i)*a_n;
end
a(n)  =s(n-1)*b(n-1)+c(n-1)*a_n;
b(n-1)=c(n-1)*b(n-1)-s(n-1)*a_n;

% Apply Givens rotations to rhs
x=rhs;
for i=1:n-1,
  x_i=x(i);
  x(i)=c(i)*x_i-s(i)*x(n);
  x(n)=s(i)*x_i+c(i)*x(n);
end

% Triangular solve
x(n)=a(n)\x(n);
x(n-1)=a(n-1)\(x(n-1)-b(n-1)*x(n));
for i=n-2:-1:1,
  x(i)=a(i)\(x(i)-d(i)*x(n)-b(i)*x(i+1));
end

% Check the answer
res_nrm_cycl_bidiag_QR=norm(M*x-rhs)
diff=norm(x_full_matrix-x)

Obviously, one should translate this Octave/Matlab code into Fortran or C etcetera, if computational efficiency is important. Even for small $n$, for example, $n=100$, such an implementation of this method maybe much faster than a general sparse $LU$ factorization or a dense LAPACK/BLAS $LU$ factorization.

$\endgroup$
0
$\begingroup$

How small is "small/moderate"? N=10? N=10000? N=1000000? How many of these systems of equations do you have to solve (one? a few hundred? a few million?)

If N is small and the total number of these systems that you have to solve is small, then there's no point in doing anything more sophisticated than using a dense LU factorization ala LAPACK/BLAS.

If N is really big and and the total number of these that you have to solve is small, then use a sparse LU factorization routine to do the job.

Only if N is really big and you have lots of these to solve would there be any point in trying to develop a custom solution to the problem.

$\endgroup$
  • $\begingroup$ $n=10$ or $n=100$ seems reasonable for our application, and I have to solve several of them (let's say $10000$) but each with a different matrix. It seems that one can do better than vanilla LU factorization; for instance, the matrix is lower Hessenberg, so one can compute its QR factorization in $O(n^2)$. And I was hoping some FFT-based trick could yield linearithmic time. $\endgroup$ – Federico Poloni Oct 21 '16 at 7:53
  • 1
    $\begingroup$ I disagree with this attitude "if it is not large scale it doesn't matter". There is a point in solving an easy small-scale problem in 0.05 seconds instead of 0.5, if it later ends up as the computational core of a longer task. $\endgroup$ – Federico Poloni Oct 21 '16 at 7:57

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.