Comments (28)
This post on stackoverflow.com might provide a starting point:
from ojalgo.
But that starting point highlights the problem of big exponents; won't it be better if we kept numerator as a pair of exponent and mantissa? It's still a proper rational number, and all the operations are well defined for it. Admittedly, we get to the land of having two kinds on rational numbers, and growing complexity, but it doesn't sound as a particularly high price to pay.
from ojalgo.
And the big-exponents problem isn't even that tricky to reproduce; if we try 1e7
as a parameter, we get the value of 5368709120000000 * (1<<23) / 4503599627370496
and the numerator is already overflown.
from ojalgo.
Dividing by gcd fixes the problem with 1e7
specifically, but leaves a question of testing:
- what are the "important" values?
- do we want
0.1
to be(3602879701896397 / 36028797018963968)
or(1 / 10)
?
from ojalgo.
I'd say the resulting rational number should be as exact as possible.
For the valueOf method 0.1 should definitely be represented as (1 / 10). Hard to explain why it should be anything else.
For the add/subtract/multiply/divide methods efficient execution is more important.
Do you have an implementations to test?
from ojalgo.
The implementation, which is quite naive, is at https://github.com/optimatika/ojAlgo/compare/master...alf239:issue49?expand=1; it fails some 233 tests (which I guess answers some of my questions), so I wasn't planning to push it just yet.
As for 0.1d
, it's 11111110111001100110011001100110011001100110011001100110011010
, which is
- sign:
0
, positive - exponent:
1111111011
(subtracting the bias of 1023, we get 0b1111111011 - 1023 = -4) - fraction:
1001100110011001100110011001100110011001100110011010
standing for binary (sic!) 1.1001100110011001100110011001100110011001100110011010 which is 0b11001100110011001100110011001100110011001100110011010 * 2**(-52) == 7205759403792794 / 4503599627370496
Now 7205759403792794 and 4503599627370496 have a GCD of 2, so while 7205759403792794 / 4503599627370496 == 1.6 in the normal world, our precise world tells us that it's actually 3602879701896397 / 2251799813685248, and if we account for the exponent (2**-4), we get 3602879701896397 / 36028797018963968, which is almost but not quite equal to 0.1:
>>> 3602879701896397.0 / 36028797018963968
0.1
I will take a look at the BigDecimal's trick around 0.1, but the naïve approach leads us to the fraction above.
from ojalgo.
Updated to use approach from https://stackoverflow.com/a/13224190/562388 — it's less readable, but follows the javadoc, which is a great thing.
from ojalgo.
Looks clean.
I don't follow the logic entirely: how many test cases fail now?
I assume to minimise the number of failing tests one should try to mimic the behaviour of the BigDecimal/BigInteger based implementation - the part with the shift.
from ojalgo.
There are 233 failing tests at the moment; I will look closer today or tomorrow.
e.g.
# CompareImplementations.testASINH
junit.framework.AssertionFailedError: Quaternion vs Rational
Expected :-0.0711857704764945
Actual :-0.07665118472290924
# CompareImplementations.testMULTIPLY
junit.framework.AssertionFailedError: Quaternion vs Rational
Expected :-0.8976069672109909
Actual :0.8976069672109909
The latter gives hope that there's some ridiculous bug hiding, stay tuned.
from ojalgo.
Regarding the number of failing tests, it's hard to tell: the tests are randomised, and the seed isn't fixed, so what fails one time can easily pass next time.
from ojalgo.
Oh, I see, so you jump to BigIntegers and then divide both numerator and denominator by 2 until they fit safely? Interesting, that may help in many cases.
from ojalgo.
Why do you use
if (retDenom.signum() > 0) {
retNumer = retNumer.negate();
retDenom = retDenom.negate();
}
shouldn't it be
if (retDenom.signum() < 0) {
retNumer = retNumer.negate();
retDenom = retDenom.negate();
}
instead?
from ojalgo.
The above was indeed a bug responsible for CompareImplementations.testMULTIPLY
failures and quite a few others. Fixed now, checking the rest.
from ojalgo.
org.ojalgo.matrix.BasicMatrixTest#testGetCondition
apparently exposes another problem: sometimes,
public Array1D<Double> getSingularValues() {
if ((mySingularValues == null) && this.isComputed()) {
mySingularValues = this.makeSingularValues();
}
return mySingularValues;
}
appears to be in situation where mySingularValues == null
and not this.isComputed()
, so the result is null
, and the test explodes with NPE. Wonder why it isn't computed...
from ojalgo.
So it's there as the matrix is not positive definite, thus cannot be decomposed.
from ojalgo.
Also, RationalNumber.POSITIVE_INFINITY.add(RationalNumber.POSITIVE_INFINITY)
throws an exception. Similar problem for NaN
, too.
from ojalgo.
Is the problem the zero denominator in POSITIVE_INFINITY and NaN? May have to review parts of the code to make sure it handles zero denominators.
The condition number is derived from the singular values decomposition. That should work regardless of "positive definite" or not.
That must be a bug related the retDenom.signum()
. The sign of the rational number should always be encoded in the numerator.
from ojalgo.
The problem is with POSITIVE_INFINITY
specifically. POSITIVE_INFINITY
+ POSITIVE_INFINITY
= 1 / 0
+ 1 / 0
= 2 / 0
= throw new ArithmeticException()
NaNs create a problem with doubleValue
which is calculated via BigDecimal
, and
The BigDecimal class provides no representation for NaN, +∞ or -∞.
(https://stackoverflow.com/a/28066001/562388)
from ojalgo.
Funnily enough, one of the failures is in the CompareImplementations.testMIN
Now the actual case it fails on is,
junit.framework.AssertionFailedError: Big vs Complex, -0.9645598407999487, 0.7898303544540544
Expected :-0.9645598407999487
Actual :0.7898303544540544
And it fails because for complex numbers we compare "norms", which are insensitive to the sign.
from ojalgo.
Raised PR #74
from ojalgo.
I'm surprised that tests pass repeatedly on the CI server, but so be it. Locally, I have three tests that fail regularly (though not always, it seems).
from ojalgo.
A better approach would be to use rational fraction aproximation, as in https://uk.mathworks.com/help/matlab/ref/rat.html?s_tid=gn_loc_drop - after all, we're not that bound to powers of 2 and 10.
Alternatively, we can keep the binary exponent around. That would let us keep the whole accuracy of double and then some.
from ojalgo.
It's a painful decision by the way. If we go with the binary exponent, we match doubles perfectly, but then 0.1 will be 7205759403792794 / 4503599627370496
just because of the limits of the mantissa.
On the other hand, if we go with the rational approximation, we'll have that 1 / 10
(haven't checked), but we'll lose any number greater than Long.MAX_VALUE
, i.e. 2^63-1.
Do we want to use both and fall back to the exact representation when we fall out of range? Do we want to follow Mathematica's route and go with BigRationals instead?
from ojalgo.
Everything in the org.ojalgo.scalar package should be primitives based. RationalNumber was previously implemented using BigInteger. We're moving away from that.
Having 0.1 represented as 7205759403792794 / 4503599627370496
seems a horrible solution. Don't you agree?
Rational approximation is really interesting (perhaps combined with some way to store magnitude separately). I think it's important to not be too clever here. When users choose to work with RationalNumber it is important that they get something that "is" a rational number.
from ojalgo.
Good point regarding using primitives, fair enough.
As for 0.1, I'm a bit puzzled here. I admit that the expectation is indeed to have 0.1 equal to 1/10 - people tend not to know that 0.1 is actually 0b1.1001100110011001100110011001100110011001100110011010 * 2^{-4}, and are not likely to appreciate us forcing this knowledge on them.
Still, every now and again that abstraction leaks, for example for 0.1 + 0.2:
>>> 0.1 + 0.2
0.30000000000000004
>>>
Mathematica deals with that problem by clearly separating the exact math (including rationals) and the approximate math; multiplying a rational by double gives one a double, not a rational.
In[1]:= 1 / 10
Out[1]= 1/10
In[2]:= 0.1
Out[2]= 0.1
In[4]:= Rationalize[0.1]
Out[4]= 1/10
In[5]:= 0.1 * (1/10)
Out[5]= 0.01
In[6]:= Rationalize[0.01]
Out[6]= 1/100
In[7]:= 0.1 * 10
Out[7]= 1.
In[8]:= (1/10) * 10
Out[8]= 1
In[9]:= 0.1 / (1/10)
Out[9]= 1.
If we follow this approach, then we should probably make the rational approximation a required call and retire the methods in the original post, as they invite the user in the land of approximation even though we promised it would be a rational.
from ojalgo.
OK, I think maybe we should try do disregard how things "look" and focus on what works best in terms of the computations.
Regarding large numbers: When the denominator gets large (compared to the numerator) then it doesn't seem problematic to approximate this with "zero". But, when it is the numerator that gets large then maybe sometimes returning the largest possible number (denominator 1) is preferable to returning infinity (denominator 0). Once you return infinity or NaN all further operations are defunct.
BTW I want to release v45 as soon as possible.
from ojalgo.
Let's see how the rational approximation approach works. I've drafted it as described, no optiization, nothing, and it seems that so far it... just works.
Once you return infinity or NaN all further operations are defunct.
But then if I don't, the result is incorrect, and heavily biased. For example, instead of 2^100 - 2^102 == -3/4 * 2^102 we would get Long.MAX_VALUE - Long.MAX_VALUE = ... well, either 0 or Long.MAX_VALUE depending on how exactly we do the subtraction.
from ojalgo.
Run complete. Total time: 00:01:21
Benchmark Mode Cnt Score Error Units
RationalCreation.rational thrpt 20 81612.897 ± 1092.180 ops/s
RationalCreation.valueOf thrpt 20 132289752.456 ± 10159005.133 ops/s
from ojalgo.
Related Issues (20)
- how to use Expression to form expression like x/1-x? HOT 1
- Solution does not balance based on weight HOT 3
- make ExpressionsBasedModel.isInfeasible() public HOT 13
- Expression should have "add" method that takes Variable as argument. HOT 4
- Eigenvalues and Eigenvectors differing from MatLab implementation? HOT 3
- getEigenpairs() doesn't work on hermitian complex matrix HOT 2
- Generalised Eigenvalue Problems do not work on complex matrices HOT 4
- Array index out of bounds exception while using IntegerStrategy HOT 7
- Non linear programming HOT 2
- hardware profile HOT 6
- Add a bibtex reference for citing this software HOT 1
- Do you have the ojAlgo js library HOT 1
- Get UNEXPLORED status after sometime, HOT 4
- Request: Support Implicit Expansion with Arithmetic Operations HOT 2
- Issue with Cholesky decomposition - symmetric positive-definite matrix decomposition is not solvable HOT 2
- invert method got wrong result HOT 3
- Solver returns INFEASIBLE on a simple problem with a clear solution HOT 6
- Allocation of number of non-zero elements in SparseStore HOT 4
- QR Decomposition's getInPlace returns null HOT 4
- ConvexSolver didn't found correct solution HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from ojalgo.