Singular value decomposition
This exercise continues with a more challenging task in numerical programming. That is, we ask you to implement a method that computes a singular value decomposition (SVD) of a given (real) matrix.
In more precise terms, given an m × n matrix A with real entries as input, you are asked to compute a decomposition A = UΣVT where
1. U is an m × m real orthogonal matrix (UUT = I where I is the m × m identity matrix),
2. Σ is an m x n nonnegative (rectangular) diagonal matrix, and
3. V is an n x n real orthogonal matrix (VVT = I where I is the n x n identity matrix).
The three-tuple (U, Σ, V) of matrices is a singular value decomposition (SVD) of the matrix A.
Remark. To simplify the task you can assume that m, n ≤ 20. Furthermore, we allow your code to be numerically quite sloppy. Namely, if A has independent and identically distributed Gaussian entries (zero mean, unit variance), your code should return with high probability an SVD (U, E, V) such that the properties (i) A = UΣVT, (ii) UUT = I, (iii) Σ is nonnegative rectangular diagonal, and (iv) VVT = I all hold within entrywise additive error of € = 106. In practice this holds if your code passes the unit
tests.
The singular value decomposition enables easy access to many basic properties of a matrix such as its inverse, pseudoinverse, null space, range, low-rank approximations, and so forth. (See here for more.)
In the context of statistics and machine learning, the singular value decomposition is used in, for example, principal component analysis.

Hints
To compute the SVD of A it is sufficient to solve eigenvector/eigenvalue problems for the real symmetric matrices A transpose(A) and transepose(A) A.
Algorithms such as the QR algorithm or the Jacobi algorithm can be used to solve for eigenvectors and eigenvalues.

Use Scala 3. Implement the function def svd: (Matrix, Matrix, Matrix) Replace the symbol "???" with your solution. Do not use external libraries like breeze.

package svd

/* A simple class for matrix arithmetic. */

class Matrix(val m: Int, val n: Int):
  require(m > 0, "The dimension m must be positive")
  require(n > 0, "The dimension n must be positive")
  protected[Matrix] val entries = new Array[Double](m * n)

  /* Convenience constructor for square matrices */
  def this(n: Int) = this(n, n)

  /* Access the elements of a matrix Z by writing Z(i,j) */
  def apply(row: Int, column: Int) =
    require(0 <= row && row < m)
    require(0 <= column && column < n)
    entries(row * n + column)
  end apply

  /* Set the elements by writing Z(i,j) = v */
  def update(row: Int, column: Int, value: Double) =
    require(0 <= row && row < m)
    require(0 <= column && column < n)
    entries(row * n + column) = value
  end update

  /* Gives a string representation */
  override def toString =
    val s = new StringBuilder()
    for row <- 0 until m do
      for column <- 0 until n do
        s ++= " %f".format(this(row,column))
      s ++= "\n"
    s.toString
  end toString

  /* Returns the transpose of this matrix */
  def transpose =
    val result = new Matrix(n,m)
    for row <- 0 until m; column <- 0 until n do
      result(column, row) = this(row, column)
    result
  end transpose   

  /* Returns a new matrix that is the sum of this and that */
  def +(that: Matrix) =
    require(m == that.m && n == that.n)
    val result = new Matrix(m,n)
    (0 until m*n).foreach(i => { 
      result.entries(i) = this.entries(i) + that.entries(i) 
    })
    result
  end +

  /* Returns a new matrix that negates the entries of this */
  def unary_- =
    val result = new Matrix(m,n)
    (0 until m*n).foreach(i => { 
      result.entries(i) = -entries(i)
    })
    result
  end unary_-

  /* Returns a new matrix that is the difference of this and that */
  def -(that: Matrix) = this + -that

  /* Returns a new matrix that is the product of 'this' and 'that' */
  def *(that: Matrix): Matrix =
    require(n == that.m)
    val thatT = that.transpose  // transpose 'that' to get better cache-locality
    def inner(r1: Int, r2: Int) =
      var s = 0.0
      var i = 0
      while i < n do
        // the inner loop -- here is where transposing 'that' pays off
        s = s + entries(r1+i)*thatT.entries(r2+i) 
        i = i+1
      s
    val result = new Matrix(m,that.n)
    (0 until m*that.n).foreach(i => { result.entries(i) = inner((i/that.n)*n, (i%that.n)*n) })
    result
  end *

  /* 
   * YOUR TASK:  
   * Implement the following method that returns the singular value
   * decomposition of this matrix.
   *
   */

  /* Returns the singular value decomposition of this matrix */
  def svd: (Matrix, Matrix, Matrix) =
    ???
  end svd

end Matrix

object Matrix:
  val r = new util.Random(1234567)

  def rand(m: Int, n: Int) =
    val a = new Matrix(m,n)
    for i <- 0 until m; j <- 0 until n do { a(i,j) = r.nextDouble() }
    a

  def randGaussian(m: Int, n: Int) =
    val a = new Matrix(m,n)
    for i <- 0 until m; j <- 0 until n do { a(i,j) = r.nextGaussian() }
    a

  def identity(n: Int) =
    val a = new Matrix(n)
    (0 until n).foreach(i => a(i,i) = 1)
    a
end Matrix