-
-
Notifications
You must be signed in to change notification settings - Fork 41
Math DHB wk
mainly adds some methods to Math-DHB-Numerical, hence it obviously depends on that package.
there are class methods to produce random vectors and matrices. you have to define a maximum number and the result will be a vector or matrix with numbers between 0 inclusively and the maximum number exclusively. if you enter an integer the result will consist of integers. examples:
DhbVector new: 7 random: 2 . "-->
a DhbVector(1 1 0 1 0 1 1)"
DhbVector new: 2 random: 2.0 ." -->
a DhbVector(1.6048663263716363 1.079899274069224)"
DhbMatrix rows: 2 columns: 3 random: 4. "-->
a DhbVector(2 1 3)
a DhbVector(3 0 0)"
DhbMatrix rows: 2 columns: 1 random:1.0 . "-->
a DhbVector(0.2250882150355371)
a DhbVector(0.49903106230149663)"
DhbSymmetricMatrix new: 3 random: 9 ."-->
a DhbVector(4 8 0)
a DhbVector(8 7 6)
a DhbVector(0 6 8)"This is defined in DHB-Numerical for Numbers and is essentially the same as closeTo:, that is it tests whether two floating point numbers are 'equal for practical purposes'. it is more rigorous than closeTo:, but as a compensation it tests the floating point precision of your computer first and adjusts its testing accordingly. btw this test is only done once, hence repeated uses of equalsTo: are not really slow. i find this method quite useful, hence i implemented it also for SequenceableCollection and DhbMatrix. examples:
#(1 3)equalsTo: #(1.00000001 3). "-->true"
#(1 3)equalsTo: #(1.00000002 3). "-->false"
#(1 3)equalsTo:#(3 1). "-->false"
a:=DhbSymmetricMatrix new: 3 random: 9.
b:=DhbMatrix rows: a rows.
a equalsTo: b. "-->true"
a:=DhbSymmetricMatrix new: 3 random: 9.0.
a equalsTo: b. "-->false"
a:=DhbSymmetricMatrix new: 2 random: 9.0.
a equalsTo: b. "--> Error: otherCollection must be the same size"QR factorization of a m-by-n Matrix A is given by
A=Q*R
with m-by-m Q being orthogonal and m-by-n R being an upper triangular matrix.
attention: the calculations in this package assume that m>=n; if this is not the case, an Error will occur!
for further info see https://en.wikipedia.org/wiki/QR_decomposition . the factorization is done via householder transformation and follows: Matrix Computations third edition by Gene H Golub & Charles F Van Loan.
DhbVector>>householder returns an Array of the skalar beta and the householder vector as described in the mentioned book.
DhbMatrix>>qrFactorization returns an Array of Q and R.
DhbMatrix>>qrFactorizationWithPivoting returns an Array of Q,R and a pivoting matrix stored in a vector as described in the book chapter 5.4.1.
for using this pivoting matrix there is the method
DhbMatrix>>inversePivotColumns:aPivotingMatrix
its use can be seen in QRTests>>testQRFactorization.
and then there are a few methods that use this factorization:
DhbMatrix>>rank calcs the rank of a matrix, something that was missing in DHB-Numerical.
DhbMatrix>>orthogonalize which returns an orthonormal basis of column vectors for a matrix of column vectors.
DhbMatrix>>mpInverse which calcs the Moore Penrose Inverse. a Moore Penrose Inverse X of matrix A satisfies these equations:
AXA=A
XAX=X
(AX)transpose=AX
(XA)transpose=XA
these 4 equations are used in QRTests>>mpTestFunction:aMatrix for testing.
DHB already has a method to calculate the Moore Penrose Inverse, but the author says in the book that this method can have significant rounding errors. additionally his method often runs into singular matrix errors. so lets do a little comparison:
a:=DhbMatrix rows: 4 columns: 3 random: 5.0.
a inverse*a. "-->
a DhbVector(0.9999999999999729 1.2434497875801753e-14 2.7373936450914016e-15)
a DhbVector(-3.019806626980426e-14 0.9999999999999454 -3.509692536596276e-14)
a DhbVector(4.218847493575595e-15 -1.4432899320127035e-14 0.9999999999999905)"
a mpInverse *a. "-->
a DhbVector(1.0000000000000009 8.881784197001252e-16 2.706168622523819e-16)
a DhbVector(8.881784197001252e-16 1.0000000000000009 5.273559366969494e-16)
a DhbVector(-4.440892098500626e-16 -2.220446049250313e-16 1.0)"
"thats obviously better"
a:=DhbMatrix rows: #((1 2 3)(3 4 6)).
a*a inverse."-->Error: ZeroDivide"
a*a mpInverse . "-->
a DhbVector(1.0 -3.3306690738754696e-16)
a DhbVector(8.881784197001252e-16 0.9999999999999993)"
"not too bad either"iterating algorithms often calculate a fixpoint, where f(x) = x by repeatedly applying f( ) on its own output. if i want to test some algo like that, i put the algo in a Block and test that block with Fixpoint. eg for the golden ratio we can repeatedly calc
xi+1 = ( 1 + xi ) / xi
lets try that:
fp:=Fixpoint block: [:x| (1 + x) / x ] value: 11.3.
"--> a Fixpoint([ :x | (1 + x) / x ]value: 11.3)"
fp evaluate .
"--> 1.618033988749895"since Fixpoint is mainly intended for simple interactive use, it returns the used number of iterations and other info via GrowlMorphs (you can shut off these GrowlMorphs by using fp verbose: false ). and of course you can also get that info 'by hand':
fp iterations.
"-->39"the default setting for the maximum number of iterations is 200, but you can reset that
fp maximumIterations: 20.
fp value: 20.0.
fp evaluate .
"-->1.6180339974620197"
"a GrowlMorph informs you that the algo has not yet converged"
fp hasConverged.
"-->false"
"but you can use #evaluate repeatedly"
fp evaluate .
"-->1.618033988749895"
fp hasConverged.
"-->true"
"btw, if you would enter an Integer as value, fp could not converge, but
calc ever more exact Fraction approximations to the irrational
golden ratio"
"lets reset things:"
fp maximumIterations: 200.
fp value: 20.0.
"and lets change the block slightly:"
fp block:[:x| (2 + x) / x ].
fp evaluate .
"-->1.9999999999999996"
"a GrowlMorph informs you that this algo has entered a 2-cycle
after 55 iterations"
fp hasConverged.
"-->true"
fp cycle.
"-->#(2.0000000000000004 1.9999999999999996)"
fp iterations.
"-->55"2-cycles are the most common problem, at least with numerical algos. for example if they switch between a Float and its successor when an irrational number is approximated. hence Fixpoint finds those cycles directly (and faster than other cycles), but it also finds cycles up to a length of maximumIterations. an example:
fp:=Fixpoint block: [:x| |d| d:=x-1. d< -2 ifTrue:[d:=0.0].d] value: 20.0.
fp evaluate .
"-->0.0"
"a GrowlMorph informs you of a 3-cycle"
fp hasConverged.
"-->true"
fp cycle.
"--> #(-1.0 -2.0 0.0)"
fp iterations.
"--> 23"of course one is not restricted to use only numerical algos, Fixpoint can use any block that can eat what he spits out:
fp block:[:x|x not ].
fp value:false.
fp evaluate .
"-->false"
"-->GrowlMorph(2-cycle, 2 iterations)"
fp hasConverged.
"-->true"
fp cycle.
"-->#(true false)"for a few examples how Fixpoint can be put to good use see eg Math-AutomaticDifferentiation.