Home > Archive > Fortran > December 2006 > one element in the inverse matrix
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
one element in the inverse matrix
|
|
| yuanliu1@gmail.com 2006-12-21, 7:05 pm |
| Hi, Everyone,
Today when I talked with my groupmate, an old topic poped up so I want
to s help here again.
I have a dense boundary element matrix A, with very large dimension.
The matrix equation is Ax=b. For a specific problem, b is all zero
except the first element. So x=A^-1_{*1}. But I only want x_1, so
basically I want A^-1_{11}.
What we are doing now is solve Ax=b by CG. but it looks redundant
because now I have all the x solved, but I only want one element. If we
solve by LUD, it looks even more clumsy because now I have all the A^-1
but I only need one element.
Are there any fast algorithm that I can solve this problem more
efficiently?
Any comments is welcome.
Rgs
Yuan
| |
| glen herrmannsfeldt 2006-12-21, 7:05 pm |
| yuanliu1@gmail.com wrote:
> I have a dense boundary element matrix A, with very large dimension.
> The matrix equation is Ax=b. For a specific problem, b is all zero
> except the first element. So x=A^-1_{*1}. But I only want x_1, so
> basically I want A^-1_{11}.
I am sure someone will tell me if I am wrong, but I believe
that (A^-1)[1,1] is A[1,1]^-1, that is, 1/A[1,1]
If you do LUD and only stop after computing the first row/column
that should also give you (A^-1)[1,1].
-- glen
| |
| Gordon Sande 2006-12-21, 7:05 pm |
| On 2006-12-21 20:28:07 -0400, glen herrmannsfeldt
<gah@seniti.ugcs.caltech.edu> said:
> yuanliu1@gmail.com wrote:
>
>
> I am sure someone will tell me if I am wrong, but I believe
> that (A^-1)[1,1] is A[1,1]^-1, that is, 1/A[1,1]
>
> If you do LUD and only stop after computing the first row/column
> that should also give you (A^-1)[1,1].
> -- glen
If the matrix is
[ a b ]
[ c d ]
then the inverse is
[ d -c ] /
[-b a ] / ( ad-bc )
so the corner element would be d / ( ad-bc ) which is not
the same as 1 / a.
All this is Cramers rule for a 2 x 2. Assuming the ASCII art works.
| |
| glen herrmannsfeldt 2006-12-21, 10:04 pm |
| glen herrmannsfeldt <gah@seniti.ugcs.caltech.edu> wrote:
> yuanliu1@gmail.com wrote:
> I am sure someone will tell me if I am wrong, but I believe
> that (A^-1)[1,1] is A[1,1]^-1, that is, 1/A[1,1]
Oops, I was a little off on that one. I was trying to remember
from doing matrix inversion by hand (for small matrices) a long
time ago.
> If you do LUD and only stop after computing the first row/column
> that should also give you (A^-1)[1,1].
OK, full inversion is O(n**3). LUD is also O(n**3), though
with a smaller coefficient. I still believe that this works, to
stop LUD after the first elements are computed, though pivoting
complicates that. LUD should be O(n) to compute each element.
If you can stop after the first one, that should be O(n) overall.
-- glen
| |
| glen herrmannsfeldt 2006-12-21, 10:04 pm |
| Gordon Sande <g.sande@worldnet.att.net> wrote:
> If the matrix is
> [ a b ]
> [ c d ]
> then the inverse is
> [ d -c ] /
> [-b a ] / ( ad-bc )
> so the corner element would be d / ( ad-bc ) which is not
> the same as 1 / a.
Yes, I figured that just after I posted the other one. Then I
had to come up with something else to say.
thanks,
-- glen
| |
| Pierre Asselin 2006-12-21, 10:04 pm |
| yuanliu1@gmail.com wrote:
> I have a dense boundary element matrix A, with very large dimension.
> The matrix equation is Ax=b. For a specific problem, b is all zero
> except the first element. So x=A^-1_{*1}. But I only want x_1, so
> basically I want A^-1_{11}.
> [ ... ]
> Are there any fast algorithm that I can solve this problem more
> efficiently?
Iterated Schur complementation ? This is like an LU factorization
where you throw away the factors as you form them. It is a bit
faster than LU factorzation but unfortunately it is still O(n^3).
Renumber your matrix so the element you need is A^{-1}_{n n}.
Assume for simplicity that your matrix is nice enough that you
can skip pivoting and get away with it. The algorithm is
do k= 1,n-1 ! pivot (k,k)
do j= k+1,n
do i= k+1,n
A(i,j)= A(i,j)-A(i,k)*A(k,j)/A(k,k)
end do
end do
end do
result= 1./A(n,n)
If you add row (or column) pivoting, you have to
exclude row n (column n) from the choice of pivot because
A(n,n) has to be the last matrix element you process.
--
pa at panix dot com
| |
| yuanliu1@gmail.com 2006-12-22, 7:05 pm |
| Hi,
Thank you very much for the suggestion. But still the gain is very
marginal. The cose of this method is
1^2+2^2+3^2+...+(n-1)^2=3D(n-1)(n)*(2n-1)/6, still O(N^3) order. It is
one step further, but not tremendous.
I would guess still a lot of redundancy.
For example, in Arnoldi method, only 1 A*x can give me the largest
eigenvalue (ideally, in reality maybe 1 is not enough). So the cost is
O(N^2), an order of smaller than O(N^3) because only very little
information is needed.
My two cents is that it is very awkward to solve the whole problem to
get just one number. I would expect O(N^2) because CG can solve the
problem with O(N^2*N_it) and get N numbers ( the solution). If I only
care about one number, I would expect some number less than that.
Are there any good way of doing it fast. The cost of Cramers rule is
evne worse than LUD. Also the matrix we are dealing with can be order
of million.
Rgds
Yuan
On 12=D4=C221=C8=D5, =CF=C2=CE=E79=CA=B121=B7=D6, p...@see.signature.invali=
d (Pierre
Asselin) wrote:
> yuanl...@gmail.com wrote:
ization[color=darkred]
> where you throw away the factors as you form them. It is a bit
> faster than LU factorzation but unfortunately it is still O(n^3).
>
> Renumber yourmatrixso theelementyou need is A^{-1}_{n n}.
> Assume for simplicity that yourmatrixis nice enough that you
> can skip pivoting and get away with it. The algorithm is
>
> do k=3D 1,n-1 ! pivot (k,k)
> do j=3D k+1,n
> do i=3D k+1,n
> A(i,j)=3D A(i,j)-A(i,k)*A(k,j)/A(k,k)
> end do
> end do
> end do
> result=3D 1./A(n,n)
>
> If you add row (or column) pivoting, you have to
> exclude row n (column n) from the choice of pivot because
> A(n,n) has to be the lastmatrixelementyou process.
>=20
> --
> pa at panix dot com
| |
| glen herrmannsfeldt 2006-12-22, 7:05 pm |
| yuanliu1@gmail.com wrote:
(snip on finding (A^-1)[1,1]
> I would guess still a lot of redundancy.
> For example, in Arnoldi method, only 1 A*x can give me the largest
> eigenvalue (ideally, in reality maybe 1 is not enough). So the cost is
> O(N^2), an order of smaller than O(N^3) because only very little
> information is needed.
> My two cents is that it is very awkward to solve the whole problem to
> get just one number. I would expect O(N^2) because CG can solve the
> problem with O(N^2*N_it) and get N numbers ( the solution). If I only
> care about one number, I would expect some number less than that.
One question, which others have asked, is if you expect to be able
to do it without pivoting. I will guess that the only way to do it
faster than O(n**3) is without pivoting, as otherwise you have to
examine too much of the matrix.
> Are there any good way of doing it fast. The cost of Cramers rule is
> evne worse than LUD. Also the matrix we are dealing with can be order
> of million.
For ones that large, you should understand the numerical problems
pretty well before trusting the answer.
-- glen
| |
| yuanliu1@gmail.com 2006-12-23, 4:42 am |
| Sorry, maybe I did not get you correctly.
In my experience of LUD, pivoting does not increase the cost much. We
can always reorder the matrix before LUD so pivoting is not a concern.
Also for the boundary value problem that I am dealing with, the
diagonal term dominates. So pivoting should not be a problem.
I know we can do LUD faster than O(N^3), but still maybe O(N^2.7) or
O(N^2.9), whatever, but the gain is marginal.
To me, it is a little bit ridiculous that I need to do almost as much
operations as LUD to get just one number.
I am not a math major so this is just an intuitive guess.
Rgds
Yuan
"glen herrmannsfeldt =D0=B4=B5=C0=A3=BA
"
> yuanliu1@gmail.com wrote:
> (snip on finding (A^-1)[1,1]
>
>
>
>
> One question, which others have asked, is if you expect to be able
> to do it without pivoting. I will guess that the only way to do it
> faster than O(n**3) is without pivoting, as otherwise you have to
> examine too much of the matrix.
>
>
> For ones that large, you should understand the numerical problems
> pretty well before trusting the answer. =20
>=20
> -- glen
| |
| Gordon Sande 2006-12-23, 10:03 pm |
| On 2006-12-23 02:06:50 -0400, yuanliu1@gmail.com said:
> Sorry, maybe I did not get you correctly.
> In my experience of LUD, pivoting does not increase the cost much. We
> can always reorder the matrix before LUD so pivoting is not a concern.
> Also for the boundary value problem that I am dealing with, the
> diagonal term dominates. So pivoting should not be a problem.
> I know we can do LUD faster than O(N^3), but still maybe O(N^2.7) or
> O(N^2.9), whatever, but the gain is marginal.
Pivotting is usually an N^2 part of an N^3 problem unless you are doing
full pivotting.
But the methods with "better" complexity often come at the price of
numerical stability. Certainly the case for Strassen's method (the 2.7
or N^log2(7) exponent.) I do not recognize the 2.9 but have seen 2.35(?)
attributed to Pan. The exponents are the asymptotes which only are effective
for truly large problems.
[color=darkred]
> To me, it is a little bit ridiculous that I need to do almost as much
> operations as LUD to get just one number.
> I am not a math major so this is just an intuitive guess.
> Rgds
> Yuan
> "glen herrmannsfeldt дµÀ£º
> "
| |
| yuanliu1@gmail.com 2006-12-23, 10:03 pm |
| Thank you very much for the comments.
Without pivoting, the proposed method above basically already solve the
whole problem, which we are currently doing.
It looks that there is no way we can solve the problem with the cost of
O(N^2).
| |
| Ron Shepard 2006-12-24, 7:06 pm |
| In article <1166907183.314569.236910@48g2000cwx.googlegroups.com>,
yuanliu1@gmail.com wrote:
> It looks that there is no way we can solve the problem with the cost of
> O(N^2).
There is no direct method that will solve the problem in O(N^2).
However, it is possible to solve linear equations using iterative
methods in O(N^2) effort. You want a single element from a single
linear equations solution, so what you s may be possible, but you
have to find an iterative method that converges quickly for your
matrix.
$.02 -Ron Shepard
| |
| yuanliu1@gmail.com 2006-12-24, 7:06 pm |
| hi, Ron
Thank you very much. Maybe direct solver is not an option.
I am using iterative solver for this problem and it converges fairly
quick, usually in less than 50 steps. Then I got the whole solution,
which are N numbers. Remebering that I only need 1 number and do not
care the others, so I expect somewhere I can speed up the process.
Iterative solver is the same as find the eigenvalues. In Arnoldi
method, one A*x theoritically should give me the largest eigenvalue,
without know the others. So if it is possible to modified the matrix so
that all the x_i are less accurate but I get very accurate evaluation
of the x_1 within very few steps?
Rgds
Yuan
| |
| Terence 2006-12-24, 7:06 pm |
| The standard method of extracting eigenvectors by finding the largest
eigenvalue first, by repeating the matrix multiplication on the trial
vector, has a potential trap. If you start with a startpoint vector
which in all coordinates is directly opposed to the true vector
"quadrant", the partial solution does not converge and oscillates about
this opposing vector. An example given in "Modern computing methods"
(HMSO 1961) actually does not work when applied.
It is the general case that all values in the matrix contribute to each
of the diagonal elements after reductions are applied. However, the
sensitivity to variations in any one value can vary immensly, so
approximations may sometimes work.
| |
| Ron Shepard 2006-12-25, 4:06 am |
| In article <1166995643.818602.55620@48g2000cwx.googlegroups.com>,
yuanliu1@gmail.com wrote:
> hi, Ron
> Thank you very much. Maybe direct solver is not an option.
> I am using iterative solver for this problem and it converges fairly
> quick, usually in less than 50 steps. Then I got the whole solution,
> which are N numbers. Remebering that I only need 1 number and do not
> care the others, so I expect somewhere I can speed up the process.
I don't know of any way to achieve what you want unless you know
more about your matrix. For example, if you know the inverse of an
(n-1)x(n-1) subblock (or, equivalently, its LU factorization), then
you can solve your problem using a direct method with only N^2
effort.
> Iterative solver is the same as find the eigenvalues.
It is roughly the same as finding a single eigenpair, not *all* of
the eigenvalues.
> In Arnoldi
> method, one A*x theoritically should give me the largest eigenvalue,
> without know the others.
Only if you know ahead of time that x is the eigenvector
corresponding to the largest eigenvalue. Otherwise, you only obtain
an upper and lower bound for *some* eigenvalue, not necessarily the
largest one.
$.02 -Ron Shepard
| |
| NuclearWizard 2006-12-25, 7:06 pm |
| > Are there any fast algorithm that I can solve this problem more
> efficiently?
>
You gave some very nice information about the rhs, namely that only the
first element is nonzero. I believe you should look for ways to
exploit this fact within an iterative method.
Write out the equations for a couple iterative methods and take into
account that only 1 term in b is nonzero and hopefully you'll see some
simplification that can trim a significant number of operations from
your solution. The downside is you have to code your own iterative
solver, but there's always a tradeoff! Good luck. :)
+Will
Example:
You may see an advantage for your special system in moving to a
"simple" iterative method and considering your explicit b. Consider
the simplest, Jacobi.
xi(k)=( bi - sum(aij*xj(k-1),j/=i) )/aii
You only care about x1, so one nice thing is that for each iteration,
because you only need the previous iterate, you can check convergence
of x1 before you calculate the rest. Also, before solving you can
modify each row i of your system by dividing by aii. Then you can
iterate on
xi(k) = bi - sum(aij*xj(k-1),j/=i),
where the pre-dividing has been taken into account. So this would be
my solution:
xnew(1) = b1 - DOT_PRODUCT( A(1,:) , xold )
IF xnew(1) SATISFIES CONVERGENCE CRITERA STOP
DO i=2,N
xnew(i) = - DOT_PRODUCT( A(i,:) , xold )
END DO
xold = xnew
Of course, there would be the actual iterations loop around the pseudo
code above. You would save some storage on rhs b, only storing the one
nonzero value. As far as operation counts, you save 2*N, which
probably wouldn't really do that much for you, but this could serve as
a decent baseline to compare how much better/worse other options fare
in terms of memory/runtime requirements.
| |
| yuanliu1@gmail.com 2006-12-27, 7:04 pm |
| Thank you for the good suggestion.
Maybe I need to take another look at the "matrix computation" to see if
there is some iterative solver can make use of this properity.
Happy new year.
Rgds
Yuan
|
|
|
|
|