Mark Clements wrote:
>
> Dear all,
>
> In writing a routine for QR decomposition of matrices (http://rosettacode.o
> rg/wiki/QR_decomposition#Axiom), I needed a domain that combined Fields (+,
> -,*,/), RadicalCategory (sqrt - for normed length of a vector) and
> OrderedRing (sign). I satisfied this using:
>
> if not(R has OrderedRing) and R has ConvertibleTo(Float) then
> sign(x:R):Integer == sign(convert(x)$R @ Float)
>
> Pleasantly, this allows a QR decomposition using Integers, Fractions or Alg=
> ebraicNumbers (as well as Float and DoubleFloat):
>
> m := matrix [[12, -51, 4], [6, 167, -68], [-4, 24, -41]];
> qr m
>
> + 6 69 58 +
> |- - --- --- |
> | 7 175 175 |
> | | +- 14 - 21 14 +
> | 3 158 6 | | |
> (2) [q= |- - - --- - ---|,r= | 0 - 175 70 |]
> | 7 175 175| | |
> | | + 0 0 - 35+
> | 2 6 33 |
> | - - -- -- |
> + 7 35 35 +
>
> Type: Record(q: Matrix(AlgebraicNumber),r: Matrix(AlgebraicNumber))
>
> Is there a better solution to this problem?
Well, sign for AlgebraicNumber is less defined than for complex
numbers:
(6) -> i2 := sqrt(-2)
+---+
(6) \|- 2
Type: AlgebraicNumber
(7) -> convert(i2)@Float
>> Error detected within library code:
cannot retract nonconstant polynomial
Actually, for QR decomposition you probably should use complex
sign, that is:
sign(x) ==
x = 0 => error "sign of 0"
x/abs(x)
(you may also define sign(0) = 0, but AFAICS in QR decomposition
you should only take nonzero pivots (0 is an error)).
Of course, now you have problem how to compute abs...
Also, in complex case you need to decide if you want orthognal
Q or unitary Q? AFAICS your code gives the first, but usually
in numerical computation people want the second (which
requires hermitian transpose instead of normal transpose).
--
Waldek Hebisch
[email protected]
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/fricas-devel?hl=en.
