0.1. Approximate and Exact Data

We start out generating random numbers to define (A,b,c).

>	n := 5:  # the number of variables

>	A := RandomMatrix(n,n,generator=-1.0..1.0,outputoptions=[shape=symmetric]):

>	exact_A := convert(A,rational):

>	b := RandomMatrix(n,1,generator=-1.0..1.0):

>	exact_b := convert(b,rational):

>	c := stats[random,uniform[-1.0,1.0]](1):

>	exact_c := convert(c,rational):

>	v := seq(x[i],i=1..n):  # sequence of variable names

>	X := Vector(n,[v]):     # vector of variable names, X is capitalized

>	p := expand(Transpose(X).A.X) + (Transpose(b).X)[1] + c;

>	exact_p := expand(Transpose(X).exact_A.X) + (Transpose(exact_b).X)[1] + exact_c;

>	representation_error := exact_p - p

0.2. Functions to Evaluate Quadrics

We can evaluate the quadric as a fully expanded polynomial, or in its matrix form.

We use a random point to test our evaluation procedures:

>	r := stats[random,uniform[-1.0,1.0]](n):

>	exact_r := op(convert([r],rational)):

>	R := Vector(n,[r]):

>	exact_R := convert(R,rational):

The little r is just a sequence of numbers, used to evaluate the quadratic form as a polynomial function (done by f1 and its exact counterpart exact_f1 below).  The big R is a vector and is input to the evaluation functions which use matrix-vector arithmetic (done by f2 and exact_f2 below).

>	f1 := unapply(p,v);

>	exact_f1 := unapply(exact_p,v);

>	f1(r); exact_f1(exact_r);

>	f2 := X -> Transpose(X).A.X + (Transpose(b).X)[1] + c;

>	exact_f2 := X -> Transpose(X).exact_A.X + (Transpose(exact_b).X)[1] + exact_c;

>	f2(R); exact_f2(exact_R);

We notice no difference in the exact numbers, as mathematically, both ways to evaluate the quadric are equivalent.

Assignment One

Give high level algorithms for the functions f1 and f2, using pseudo-code making the additions and multiplications explicit.

Count the number of arithmetical operations it takes to evaluate f1 and f2.  Justify your count referring to the algorithmic descriptions of f1 and f2. Start with n=3, and continue with n=4, n = 5, etc,  till you find a general formula in terms of n.

Which one is most efficient, f1 or f2?

2.1. Evaluation at Random Points

A random point is generated first as a sequence of random numbers (input to f1), then converted into a sequence rational numbers (input to exact_f1) and into a vector of numbers (input to f2).

>	r := stats[random,uniform[-1.0,1.0]](n);

>	exact_r := op(convert([r],rational));

>	R := Vector(n,[r]);

>	y1 := f1(r); y2 := f2(R); y0 := exact_f1(exact_r);

To compute the relative errors, we temporarily increase the working precision:

>	Digits := 16;

>	relative_error1 := abs(y1-y0)/abs(y0); relative_error2 := abs(y2-y0)/abs(y0);

>	Digits := 8;

2.2. Evaluation at Points of Large Magnitude

In this part we keep the n fixed to 5, but let the magnitudes of the points decrease or increase.  For example, to get a random point of magnitude 10^4, we compute

>	r4 := 10^4*r: exact_r4 := op(convert([r4],rational)): R4 := Vector(n,[r4]):

>	y1 := f1(r4); y2 := f2(R4); y0 := exact_f1(exact_r4);

Also here we raise the precision to compute the errors:

>	Digits := 16:

>	absolute_error1 := abs(y1-y0); absolute_error2 := abs(y2-y0);

>	relative_error1 := absolute_error1/abs(y0); relative_error2 := absolute_error2/abs(y0);

>	Digits := 8:

2.3. Evaluation at a Special Point

The special point of interest is the so-called center of the quadric, defined as the solution to Ax = -b.  Evaluating at the center should give us c back.

We first calculate the center exactly before rounding it to our working precision.

>	exact_S := Column(LinearSolve(exact_A,-exact_b),1):

>	exact_s := op(convert(exact_S,list)):

>	s := evalf(exact_s):

>	S := Vector(n,[s]):

>

f1(s);

>	exact_f1(exact_s) = exact_c;

>

evalf(%);

>

f2(S);

>

c;

Assignment Two -- conclusion

Based on all experiments, write a conclusion in one paragraph.