matrixfunctions / GraphMatFun.jl

Computation graphs for matrix functions
MIT License
13 stars 0 forks source link

Normalization of first degopt row #21

Closed jarlebring closed 3 years ago

jarlebring commented 3 years ago

Without loss of generality it seems the first row in the degopt form can always be normalized to (0 1)(0 1) such that it corresponds to x^2. Justification: We can merge coefficients. If first row is (a b)(c d) <=> ac+(ad+bc)*x+bd*x^2 We can multiply this out in an arbitrary row later:

(q1 q2 q3 q4 q5)(p1 p2 p3 p4 p5)
<=> 
(q1+q2*x+q3* (ac+(ad+bc)*x+bd*x^2) + ... )  x 
     (p1+p2*x+p3* (ac+(ad+bc)*x+bd*x^2) + ... )
= 
(q1+q3*ac+(q2+q3*(ad+bc))*x+q3*bd*x^2 + ... )  x 
     (p1+p3*ac+(p2+p3*(ad+bc))*x+p3*bd*x^2 + ... ) 

Therefore if we change the first row from (a b)(c d) to (0 1)(0 1) we can always take this change into account in subsequent row by:

adding to first element:q3*ac adding to second element: q3*(ad+bc) setting the third element to: q3*bd

Proposed code: We implement normalize_degopt!(graph,:row1) which transforms the graph into an equivalent code. (The :row1 can be useful since I think there is a generalization of this reasoning.)

Why?: a) It reduces the number of lincomb b) The modern implementation of the matrix exponential (by Awad et al) makes an estimate of the spectral radius / norm by computing powers which anyway have to be computed later. By normalizing the degopt we would allow for the same trick at least for the first power.

jarlebring commented 3 years ago

With a similar technique it seems the second row can be normalized to (0 0 1) (x x x) if the first row is normalized.

I also think (0 0 1) (0 0 1) should be possible, since it seems possible to compensate for x^3 terms, but the details become complicated.

jarlebring commented 3 years ago

Next step: Suppose the first row is normalized and the second row is (0 0 1)(a b c) <=> a*x^2+b*x^3+c*x^4=a*x^2+x^2*(b*x+c*x) Since the term ax^2 can be incorporated into any later use, we can normalize like this (0 0 1)(0 b c) Assuming b!=0, we can even set the b-coeff to 1. `(0 0 1)(0 1 z) <=> x^2(x+z*x^2)`

jarlebring commented 3 years ago

Row 1 implemented. Once we figure out row2 or arbitrary, we can open a new issue.