[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Full matrix inversion for triangular and
From: |
John W. Eaton |
Subject: |
Full matrix inversion for triangular and |
Date: |
Wed, 6 Dec 2006 13:11:01 -0500 |
On 6-Dec-2006, David Bateman wrote:
| The attached patch adds probing of the matrix type and special casing of
| positive definite and triangular matrices to the Matrix and
| ComplexMatrix classes inverse methods and adjusts the Finv and xpow
| functions appropriately. As an example consider
|
| N=1000;
| A = 10*eye(N) + randn(N,N);
| t0 = cputime(); BA = inv(A); t(1) = cputime() - t0; ## Normal case
| P = A*A';
| t0 = cputime(); BP = inv(P); t(2) = cputime() - t0; ## PD case
| U = triu(A);
| t0 = cputime(); BU = inv(U); t(3) = cputime() - t0; ## Upper Triangular
| L = tril(A);
| t0 = cputime(); BL = inv(L); t(4) = cputime() - t0; ## Lower Triangular
| t
|
| which returns
|
| 7.7648 7.7618 3.8974 4.0144
|
| for 2.9.9, and for 2.9.9+ with the attached patch returns
|
| 7.8098 3.7634 1.3958 1.4018
|
| for significant gains in speed. Note that this essentially makes cholinv
| redundant, though it appears that the same is not so for chol2inv. Consider
|
| t0 = cputime(); BP = cholinv(P); cputime() - t0
| ans = 3.7494
|
| R = chol(P); t0 = cputime(); BP = chol2inv(R); cputime() - t0
| ans = 2.6386
| t0 = cputime(); BP = inv(R); BP = BP * BP'; cputime() - t0
| ans = 5.3012
|
| Now to check correctness, I did
|
| N = 100;
| n = 10;
| for i=1:n
| A = 10*eye(N) + randn(N,N);
| P = A*A';
| assert(A*inv(A), eye(N), 1e-10);
| endfor
|
| for i=1:n
| A = 10*eye(N) + randn(N,N);
| U = triu(A);
| assert(U*inv(U), eye(N), 1e-10);
| endfor
|
| for i=1:n
| A = 10*eye(N) + randn(N,N);
| L = tril(A);
| assert(L*inv(L), eye(N), 1e-10);
| endfor
|
| which ran through without assert signaling any failures. This patch
| excludes the toeplitz solver code I previously sent
OK, please check in these changes.
Thanks!
jwe