|
»
»
Octave hangup for expm with non-finite arguments
|
View:
New views
10 Messages
—
Rating Filter:
Alert me
|
 |
Octave hangup for expm with non-finite arguments
by Rolf Fabian
::
Rate this Message:
octave-3.0.0.exe:35> computer, version
i686-pc-msdosmsvc ans = 3.0.0 octave-3.0.0.exe:36> -Inf*0 ans = NaN # OK by definition octave-3.0.0.exe:37> exp(-Inf) ans = 0 # OK octave-3.0.0.exe:38> expm(-Inf) ans = 0 # OK octave-3.0.0.exe:39> -Inf*eye(2) # OK ans = -Inf NaN NaN -Inf octave-3.0.0.exe:40> expm(x= -Inf*eye(2)) # FAIL HANGUP : Had to press Ctrl-C to get back Octave prompt In order to check whether this is due to the 'NaN' entries of x I replaced them by zeros, but Octave still hang up : octave-3.0.0.exe:40> expm(y=[-Inf,0;0,-Inf]) # FAIL HANGUP : Had to press Ctrl-C to get back Octave prompt In order to check what we can expect, y might be approximated by using realmax. octave-3.0.0.exe:41> z = -realmax*eye(2) x = -1.7977e+308 -0.0000e+000 -0.0000e+000 -1.7977e+308 octave-3.0.0.exe:42> expm(z) ans = 1 0 0 1 Hmm ... expm ( - Inf ) = 0 and expm ( [-Inf,0;0,-Inf] ) = eye(2) ? Is this consistent ? Have to think about this point. Of course the much more important issue concerns the hangup of Octave at line 40. Rolf Fabian
<r dot fabian at jacobs-university dot de> |
|
|
Re: Octave hangup for expm with non-finite arguments
by Rolf Fabian
::
Rate this Message:
I've analyzed this issue in the meantime. expm ( [-Inf,0;0,-Inf] ) should return zeros(2), and not ones(2) ! More generally for arguments X of expm with repeated -Inf on the main diagonal ( and zeros elsewhere ) X = diag( -Inf*ones(1,n) ) =========================================== expm( X ) should return zeros(n) for all dimensions n>0 . =========================================== Thus X can be considered to be the matrix logarithm of zeros(n) matrices. ============================== logm( zeros(n) ) = diag( -Inf*ones(1,n) ) ============================== Proof: This definition is reasonable, because it also results from the limiting process outlined below. I figured out the following general argument scaling formula for the matrix logarithm function ( non-zero real or complex scaling factor s required ): ================================== logm ( s*x ) = logm ( x ) + log (s) .* eye (n) ================================== ( it may be checked with any arbitrary non-zero x but keep in mind that Octave's logm(x) may fail if your selection of x is 'defective'. E.g. for random matrices this is extremely unlikely. ) If we apply the general scaling formula to the particular matrix x = eye(n) and look at the limiting process s -> 0, we'll get lim s-> 0 logm ( eye(n)*s ) = logm ( eye(n) ) + log (s) .* eye(n) = logm ( zeros(n) ) = zeros(n) + log(0).*eye(n) ) = = -Inf .* diag( ones(1,n) ) q.e.d Alignment with ML ? Can anybody with access to MatLab please check what it returns for the following code : -Inf * eye ( 2 ) expm( [ -Inf, 0; 0, -Inf ] ) logm ( zeros (2 ) ) Thanks in advance. Rolf Fabian < r dot fabian at jacobs-university dot de> Rolf Fabian
<r dot fabian at jacobs-university dot de> |
|
|
|
|
|
Re: Octave hangup for expm with non-finite arguments
by Ben Abbott
::
Rate this Message:
On Jan 17, 2008, at 7:28 AM, Rolf Fabian wrote:
> > Rolf Fabian wrote: >> >> octave-3.0.0.exe:40> expm(y=[-Inf,0;0,-Inf]) # FAIL >> HANGUP : Had to press Ctrl-C to get back Octave prompt >> >> In order to check what we can expect, >> y might be approximated by using realmax. >> >> octave-3.0.0.exe:41> z = -realmax*eye(2) >> x = >> -1.7977e+308 -0.0000e+000 >> -0.0000e+000 -1.7977e+308 >> >> octave-3.0.0.exe:42> expm(z) >> ans = >> 1 0 >> 0 1 >> >> Hmm ... >> expm ( - Inf ) = 0 and >> expm ( [-Inf,0;0,-Inf] ) = eye(2) ? >> Is this consistent ? >> Have to think about this point. >> Matlab gives >> expm([-Inf,0;0,-Inf]) Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. > In expm>PadeApproximantOfDegree at 123 In expm at 39 ans = NaN NaN NaN NaN >> > I've analyzed this issue in the meantime. > > expm ( [-Inf,0;0,-Inf] ) should return zeros(2), and not ones(2) ! > > More generally for arguments X of expm with repeated -Inf > on the main diagonal ( and zeros elsewhere ) > X = diag( -Inf*ones(1,n) ) > =========================================== > expm( X ) should return zeros(n) for all dimensions n>0 . > =========================================== > > Thus X can be considered to be the matrix logarithm of > zeros(n) matrices. > ============================== > logm( zeros(n) ) = diag( -Inf*ones(1,n) ) > ============================== > > Proof: > This definition is reasonable, because it also > results from the limiting process outlined > below. > > I figured out the following general argument scaling > formula for the matrix logarithm function > ( non-zero real or complex scaling factor s required ): > > ================================== > logm ( s*x ) = logm ( x ) + log (s) .* eye (n) > ================================== > ( it may be checked with any arbitrary non-zero x > but keep in mind that Octave's logm(x) may fail > if your selection of x is 'defective'. E.g. for > random matrices this is extremely unlikely. ) > > If we apply the general scaling formula to the particular > matrix x = eye(n) and look at the limiting process > s -> 0, we'll get > > lim s-> 0 logm ( eye(n)*s ) = logm ( eye(n) ) + log (s) .* eye(n) > = logm ( zeros(n) ) = zeros(n) + log(0).*eye(n) ) = > = -Inf .* diag( ones(1,n) ) q.e.d > > Alignment with ML ? > Can anybody with access to MatLab please check what > it returns for the following code : > > -Inf * eye ( 2 ) > expm( [ -Inf, 0; 0, -Inf ] ) > logm ( zeros (2 ) ) >> logm(zeros(2)) Warning: Principal matrix logarithm is not defined for A with nonpositive real eigenvalues. A non-principal matrix logarithm is returned. > In funm at 156 In logm at 27 ans = NaN NaN NaN NaN > Thanks in advance. > Regarding logm(0) being -Inf, as that unique? What I'm asking is, although exp(-Inf) = 0, is it proper to conclude that log(0) must be -Inf? It would appear to me that there is no manner to determine a unique solution. All we can be sure of is that the real part must be negative and infinite, but what of the imaginary part? ... it can take any finite or infinite value correct? I'm not familiar with how the IEEE defines such things, but it would appear to me to be a NaN, but what does the IEEE say? Ben _______________________________________________ Bug-octave mailing list Bug-octave@... https://www.cae.wisc.edu/mailman/listinfo/bug-octave |
|
|
Re: Octave hangup for expm with non-finite arguments
by Rolf Fabian
::
Rate this Message:
For scalar argument x, logm(x) and log(x) are expected to be identical. No ,matter whether x is scalar or a square matrix, the logarithm is always an function with infinite values. The whole set of values of y= log (x=0) is obtained by adding 'phase arguments' y = -Inf + k .* 2 .* pi .* i where k is an arbitrary integer (negative, positive or zero) Octave returns only principal value ( using k=0 here ) octave-3.0.0.exe:17> log(0) |- ans = -Inf exp ( y ) reconstructs argument of log, i.e. x no matter what the phase is. This can be easily seem by: octave-3.0.0.exe:20> exp( log(0) + (randperm(5)-3).*2*pi*i ) ans = 0 0 0 0 0 For square matrix logarithm (logm) the situation is similar to the scalar log - fuction. It has also infinite values. The major difference is that the phase arguments need to be square matrices of same size as input x and they must conform to matrix logarithms of idntity matrix: phase = logm ( eye (dimension n) ). The trivial principal logarithm is of course phase=zeros(n) but infinitely more do exist. So if y = logm ( x ), then with any y = y .+ (phase matrix) the call 'expm ( y )' reconstructs x. We do not need to get the whole set of solutions for log(m)-functions. The particular one which is called pricipal logarithm is enough. BTW MatLab uses the same concept. If you look at the ML error message above it says "principal logarithm not defined ..." My interpretation of this message is that using trivial (see above) phase matrix is not possible in this particular case Your statement about 'arbitrary finite and infinite' imaginary parts of a matrix logarithm is a misconclusion. My suggestion was to define a principal matrix logarithm of zeros(n) by y = diag ( -Inf * ones (1,n) ) i.e. Impaginary part should be zero. I'm not an expert of IEEE either. But I'm afraid that they do not consider such issues. BTW Do you also see that your Octave system hangs for command 'expm([-Inf,0;0,-Inf]) ? Rolf Fabian < r dot fabian at jacobs-university dot de> Rolf Fabian
<r dot fabian at jacobs-university dot de> |
|
|
Re: Octave hangup for expm with non-finite arguments
by Sergei Steshenko-2
::
Rate this Message:
--- Rolf Fabian <Rolf.Fabian@...> wrote: > > BTW Do you also see that your Octave system > hangs for command 'expm([-Inf,0;0,-Inf]) ? > I do not see it hanging (octave-3.0.0, Limux 32 bits): " octave:20> expm([-Inf,0;0,-Inf]) error: matrix singular to machine precision warning: dgelsd: rank deficient 2x2 matrix, rank = 0 ". --Sergei. Applications From Scratch: http://appsfromscratch.berlios.de/ ____________________________________________________________________________________ Never miss a thing. Make Yahoo your home page. http://www.yahoo.com/r/hs _______________________________________________ Bug-octave mailing list Bug-octave@... https://www.cae.wisc.edu/mailman/listinfo/bug-octave |
|
|
Re: Octave hangup for expm with non-finite arguments
by John W. Eaton
::
Rate this Message:
On 17-Jan-2008, Marco Caliari wrote:
| The attached patches should addres the expm(-realmax) bug. | | * Octave 3.0.0 with the patches | | octave:6> expm([-Inf,0;0 -Inf]) | error: matrix singular to machine precision, rcond = nan Should we error here, or give a warning? I fixed it to give a warning and return a matrix of NaNs. That might not be the best thing, but I don't know how to do any better here. Someone who is an expert will have to provide a patch. | There was an overflow in the scale_factor. OK, instead of repeatedly dividing the matrix m by two, I checked in the follwoing change that limits sqpow to 1023. With this patch, I see: octave:1> expm ([-Inf, 0; 0, -Inf]) warning: singular matrix encountered in expm calculation, rcond = 0 ans = NaN NaN NaN NaN octave:2> expm ([-realmax, 0; 0, -realmax]) ans = 0 0 0 0 However, there is still a problem. A calculation like expm (-realmax * eye (1000)); will take an extremely long time to complete because of the loop: // Reverse preconditioning step 3: repeated squaring. while (sqpow) { retval = retval * retval; sqpow--; } in which sqpow will be 1023. Multiplying a large matrix by itself this many times is going to be trouble. So we could use a better way of doing this. I don't have time to work on a fix for this problem, but would consider a patch. Thanks, jwe 2008-01-18 John W. Eaton <jwe@...> * dMatrix.cc (solve_singularity_warning): New function. (Matrix::expm): Pass pointer to solve_singularity_warning to Matrix::solve method. Exit early if Matrix::solve fails. Limit sqpow value to avoid overflowing scale factor. * CMatrix.cc (solve_singularity_warning): New function. (ComplexMatrix::expm): Pass pointer to solve_singularity_warning to ComplexMatrix::solve method. Exit early if ComplexMatrix::solve fails. Limit sqpow value to avoid overflowing scale factor. From Marco Caliari <mcaliari@...>. Index: liboctave/CMatrix.cc =================================================================== RCS file: /cvs/octave/liboctave/CMatrix.cc,v retrieving revision 1.143 diff -u -u -r1.143 CMatrix.cc --- liboctave/CMatrix.cc 6 Dec 2007 19:16:47 -0000 1.143 +++ liboctave/CMatrix.cc 18 Jan 2008 08:13:13 -0000 @@ -2750,6 +2750,14 @@ 1.9270852604185938e-9, }; +static void +solve_singularity_warning (double rcond) +{ + (*current_liboctave_warning_handler) + ("singular matrix encountered in expm calculation, rcond = %g", + rcond); +} + ComplexMatrix ComplexMatrix::expm (void) const { @@ -2843,6 +2851,9 @@ if (sqpow > 0) { + if (sqpow > 1023) + sqpow = 1023; + double scale_factor = 1.0; for (octave_idx_type i = 0; i < sqpow; i++) scale_factor *= 2.0; @@ -2889,8 +2900,12 @@ // Compute pade approximation = inverse (dpp) * npp. - retval = dpp.solve (npp); - + double rcond; + retval = dpp.solve (npp, info, rcond, solve_singularity_warning); + + if (info < 0) + return retval; + // Reverse preconditioning step 3: repeated squaring. while (sqpow) Index: liboctave/dMatrix.cc =================================================================== RCS file: /cvs/octave/liboctave/dMatrix.cc,v retrieving revision 1.147 diff -u -u -r1.147 dMatrix.cc --- liboctave/dMatrix.cc 14 Jan 2008 08:58:02 -0000 1.147 +++ liboctave/dMatrix.cc 18 Jan 2008 08:13:13 -0000 @@ -2378,6 +2378,14 @@ 1.9270852604185938e-9, }; +static void +solve_singularity_warning (double rcond) +{ + (*current_liboctave_warning_handler) + ("singular matrix encountered in expm calculation, rcond = %g", + rcond); +} + Matrix Matrix::expm (void) const { @@ -2463,10 +2471,13 @@ if (sqpow > 0) { + if (sqpow > 1023) + sqpow = 1023; + double scale_factor = 1.0; for (octave_idx_type i = 0; i < sqpow; i++) scale_factor *= 2.0; - + m = m / scale_factor; } @@ -2509,8 +2520,12 @@ // Compute pade approximation = inverse (dpp) * npp. - retval = dpp.solve (npp, info); - + double rcond; + retval = dpp.solve (npp, info, rcond, solve_singularity_warning); + + if (info < 0) + return retval; + // Reverse preconditioning step 3: repeated squaring. while (sqpow) _______________________________________________ Bug-octave mailing list Bug-octave@... https://www.cae.wisc.edu/mailman/listinfo/bug-octave |
|
|
Re: Octave hangup for expm with non-finite arguments
by Marco Caliari-2
::
Rate this Message:
On Fri, 18 Jan 2008, John W. Eaton wrote:
> On 17-Jan-2008, Marco Caliari wrote: > > | The attached patches should addres the expm(-realmax) bug. > | > | * Octave 3.0.0 with the patches > | > | octave:6> expm([-Inf,0;0 -Inf]) > | error: matrix singular to machine precision, rcond = nan > > Should we error here, or give a warning? I fixed it to give a warning > and return a matrix of NaNs. That might not be the best thing, but I > don't know how to do any better here. Someone who is an expert will > have to provide a patch. The result should be [0 0;0 0] (similarly to exp(-Inf)=0). A very dirty trick could be to replace -Inf with -realmax, and give a warning. > | There was an overflow in the scale_factor. > > OK, instead of repeatedly dividing the matrix m by two, I checked in > the follwoing change that limits sqpow to 1023. With this patch, I see: > > octave:1> expm ([-Inf, 0; 0, -Inf]) warning: singular matrix encountered in expm calculation, rcond = 0 > ans = > > NaN NaN > NaN NaN > > octave:2> expm ([-realmax, 0; 0, -realmax]) > ans = > > 0 0 > 0 0 I was also thinking the same. But, for your example, sqpow would be 1025. If you truncate sqpow to 1023, the Pade' approximation of the reduced matrix is not as accurate as wanted. You should instead consider something like if (sqpow > 0) { double scale_factor = 1.0; for (octave_idx_type i = 0; i < 1023; i++) scale_factor *= 2.0; m = m / scale_factor; for (octave_idx_type i = 1023; i < sqpow; i++) m = m / 2.0; } > However, there is still a problem. A calculation like > > expm (-realmax * eye (1000)); > > will take an extremely long time to complete because of the loop: > > // Reverse preconditioning step 3: repeated squaring. > > while (sqpow) > { > retval = retval * retval; > sqpow--; > } > > in which sqpow will be 1023. Multiplying a large matrix by itself > this many times is going to be trouble. So we could use a better way > of doing this. I don't have time to work on a fix for this problem, > but would consider a patch. No idea. Matlab requires a very long time, too. Marco P.S. Please, update my email address in the ChangeLog to marco.caliari@... _______________________________________________ Bug-octave mailing list Bug-octave@... https://www.cae.wisc.edu/mailman/listinfo/bug-octave |
|
|
Re: Octave hangup for expm with non-finite arguments
by John W. Eaton
::
Rate this Message:
On 18-Jan-2008, Marco Caliari wrote:
| On Fri, 18 Jan 2008, John W. Eaton wrote: | | > On 17-Jan-2008, Marco Caliari wrote: | > | > | The attached patches should addres the expm(-realmax) bug. | > | | > | * Octave 3.0.0 with the patches | > | | > | octave:6> expm([-Inf,0;0 -Inf]) | > | error: matrix singular to machine precision, rcond = nan | > | > Should we error here, or give a warning? I fixed it to give a warning | > and return a matrix of NaNs. That might not be the best thing, but I | > don't know how to do any better here. Someone who is an expert will | > have to provide a patch. | | The result should be [0 0;0 0] (similarly to exp(-Inf)=0). A very dirty | trick could be to replace -Inf with -realmax, and give a warning. It seems there should be a better fix. | I was also thinking the same. But, for your example, sqpow would be 1025. | If you truncate sqpow to 1023, the Pade' approximation of the reduced | matrix is not as accurate as wanted. You should instead consider something | like | | if (sqpow > 0) | { | double scale_factor = 1.0; | for (octave_idx_type i = 0; i < 1023; i++) | scale_factor *= 2.0; | | m = m / scale_factor; | for (octave_idx_type i = 1023; i < sqpow; i++) | m = m / 2.0; | } Given the way sqpow is calculated, I think the effective limit on it is 1025 anyway. Does the extra factor of 4 really matter? | P.S. Please, update my email address in the ChangeLog to | marco.caliari@... OK. jwe _______________________________________________ Bug-octave mailing list Bug-octave@... https://www.cae.wisc.edu/mailman/listinfo/bug-octave |
|
|
Re: Octave hangup for expm with non-finite arguments
by Marco Caliari-2
::
Rate this Message:
On Fri, 18 Jan 2008, John W. Eaton wrote:
> | If you truncate sqpow to 1023, the Pade' approximation of the reduced > | matrix is not as accurate as wanted. You should instead consider something > | like > | > | if (sqpow > 0) > | { > | double scale_factor = 1.0; > | for (octave_idx_type i = 0; i < 1023; i++) > | scale_factor *= 2.0; > | > | m = m / scale_factor; > | for (octave_idx_type i = 1023; i < sqpow; i++) > | m = m / 2.0; > | } > > Given the way sqpow is calculated, I think the effective limit on it > is 1025 anyway. Does the extra factor of 4 really matter? Probably not. Marco _______________________________________________ Bug-octave mailing list Bug-octave@... https://www.cae.wisc.edu/mailman/listinfo/bug-octave |
| Free embeddable forum powered by Nabble | Forum Help |
