This is a special homework. You have to produce a Project Report.

While discussion among students is encouraged, code writing and submission are individual.
The student’s candidate number must appear clearly on top of the project and contents must
not be copied.
Submission is on Canvas.
Electronic submission
One single PDF file. (Any other format than PDF is not accepted). No special typesetting
system is required, and you are welcome to handwrite the theoretical part, although LATEX(or
variants thereof) is recommended.
All computer print outs code or output must be clearly labelled and described. Long listings of
code and output should be avoided. A graph plot or a table is worth a thousand uncommented
numbers.
Remarks on marking
The project is marked in a “qualitative way”; ideas and critical discussion is more important
than technical prowess (althought the latter doesn’t hurt if you can show it).
You do not have to solve the full project to get full credit (try to answer 4 or 5 items). You can
focus on either analysis or coding (or both if you like).
4.1 The matrix square root
Denote by Rd ×d
the set of real-coefficient matrices. A matrix A ∈ Rd ×d
is called positive definite if and only if for each v ∈ Rd¬0
(C.1) v
∗Av > 0.
Denote by PD(Rd
) the set of all positive definite matrices. Denote by Spec A the spectrum (the
set of all eigenvalues) of A
(C.2) Spec A :=

λ ∈ C : Aw = λw for some w ∈ C
d¬0

.
Given Y ∈ PD(R2
), we consider here the nonlinear equation of finding X ∈ PD(R2
) such that
(C.3) X
2 = Y .
(a) Prove that if A ∈ PD(Rd
) then
(i) A is invertible,
(ii) the real part of any eigenvalue λ ∈ Spec A is strictly positive; i.e., reλ > 0.
(b) Consider a matrix A ∈ PD(R2
). Prove that the following map
(C.4) LA : R2×2 → R2×2
Z 7→ AZ + Z A
is linear and invertible in Z .
Hint. You may find it useful to write the Jordan canonical form of A
(C.5) A = W J W −1
where W is an (ordered) generalised eigenbasis of A and J is of one of the following
forms
(C.6) 
λ1
λ2

or 
λ 1
λ

1
with λ1 and λ2 possibly equal or not in Spec A, or λ the sole (defective) eigenvalue of
Spec A.
Alternatively, you can use the existence theorem for the solutions of Sylvester’s equation
(see, e.g., (wikipedia)) which says that the linear equation
(C.7) AZ + Z B = C ,
is uniquely solvable for Z if and only if Spec A ∩Spec(−B) = ∅.
(c) Prove the uniqueness of X in PD(R2
) satisfying(C.3).
Hint. Prove first that if X , Xe are in PD(Rd
)then also P := X +Xe ∈ PD(Rd
) and the identity
(C.8) X
2 − Xe

2

1
2
X − Xe

X + Xe


+

X + Xe

X − Xe


.
You may then want to recall (b) for LX +Xe .
(d) Write down the Newton–Raphson iteration for(C.3) in the form
(C.9) X (k + 1) =G (X (k)) for k ≥ 0
by finding the appropriate expression for G (X ) in terms of X , Y and the linear operator
LX introduced above. Explain why it is reasonable to expect convergence assuming that
the iterates are larger than the exact solution (in the sense of positive definite matrices)
X (k) − X ∈ PD(R2
).
You may want to attempt a proof when Y is symmetric (i.e., self-adjoint).
A good initial guess satisfying this is X (0) = I + Y . Can you explain this?
(e) Use MATLAB/Octave’s fsolve to solve the equation X
2 = Y . Explain how you are
using it and all your steps.
Try your approach on various examples, e.g.,
(C.10) 
1 −1
1 1
,

4 1
0 4
,

2 −1
−1 2
,


4 3 1
0 3 −1
1 1 5


but you can make up your own (and test for definite positivity using the MATLAB/Octave
function eig , or its numpy equivalent).
(f) Now try to avoid using fsolve and implement it directly. A first attempt could be to try
to solve, at the k-th iteration, that X k having being computed the linear problem for S
(C.11) X k S + S X k = X
2
k − Y
and then updating
(C.12) X k+1
:= X k + S .
In principle you could implement a linear solver for (C.11) by using the Kronecker product
and the “noodle”, but you are allowed here to use the MATLAB/Octave built in sylvester .
Once the code works, test it, and then look at the convergence rates.
(g) The following is an inexact Newton idea, credited to Visser (1937) by Higham (2008),
that replaces the Sylvester solve in (C.11) with the following:
(C.13) D S + S D = X
2
k − Y
where D =
1
2α
I is a matrix that is constant in k, for some scalar α; the art is in choosing
the right α.
Implement it, test it and calculate the experimental convergence rates, like above. Experiment with various values of α related to the matrix Y .
(h) Use the Visser idea, to show that the solution of the square root problem has a unique
postive definite solution when Y is positiv
References
[1] Nicholas J. Higham. Functions of Matrices. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008. ISBN: 978-0-89871-646-7. DOI: 10.1137/1.9780898717778.
URL: https://search.worldcat.org/en/title/185095770 (cit. on p. 2).