Matrix inversion and applications
This exercise practices your skills in numerical programming and numerical algorithms. Namely, we ask you to implement a method that computes the inverse of a matrix. For example, you may use Gaussian elimination. Alternatively, you may want to familiarize yourself with matrix decompositions such as the LU decomposition or the QR decomposition, both of which enable easy calculation of the inverse matrix.

Remark. You can assume that the matrix is invertible and well-behaved so that numerical accuracy of floating-point arithmetic does not become an issue. (However, in case you are interested, start here for more about numerical accuracy and ill-conditioned matrices.)

Matrix inversion and/or related matrix decompositions are required as subroutines in a myriad of applications.

One classical application in the context of statistics and machine learning is least squares estimation which can be used, for example, in the context of linear regression and polynomial regression. Once you have completed your matrix inversion method, we ask you to implement a method that fits a polynomial of degree d to n data points so that the sum of squares of errors on the data points are minimized.

Hints
Gaussian elimination can be found in pseudocode on the web and in essentially any numerical algorithms textbook. The same holds for the LU and QR decompositions.
Make use of the class Matrix and your matrix inverter when implementing least-squares polynomial fitting.

Use Scala 3. Complete task 2. Replace your solution with the "???" symbol in task 2.

/* 
 * This is 'matrix.scala'. 
 * See "Task 1" below
 *
 */

package matrixInverse

/* 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) : Double =
    require(0 <= row && row < m)
    require(0 <= column && column < n)
    entries(row * n + column)

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

  /* Gives a pretty-printed string representation for small matrices */
  override def toString : String =
    val s = new StringBuilder()
    if m <= 6 && n <= 6 then
      for row <- 0 until m do
        for column <- 0 until n do
          s ++= " %f".format(this(row,column))
        s ++= "\n"
    else
      s ++= "[%d-by-%d matrix -- output suppressed]".format(m,n)
    s.toString

  /* Returns the transpose of this matrix as a new matrix.*/
  def transpose : Matrix =
    val result = new Matrix(n,m)
    for row <- 0 until m; column <- 0 until n do
      result(column, row) = this(row, column)
    result
    
  /* Returns a new matrix that is the sum of this and that */
  def +(that: Matrix) : 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

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

  /* Returns a new matrix that is the difference of this and that */
  def -(that: Matrix) : 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

  /*
   * Task 1: 
   *
   * Implement the following method that returns the multiplicative 
   * inverse of this matrix. For example, you can use Gaussian elimination.
   *
   * Remark:
   * You may assume the matrix to be invertible and numerically well-behaved.
   *
   */
  // Returns a new matrix that is the multiplicative inverse of this matrix.
  def inverse: Matrix =
    require(n == m)
    // create an augmented matrix [A | I]
    val augmented = new Matrix(m, 2 * n)
    (0 until m).foreach(i => {
      (0 until n).foreach(j => {
        augmented(i, j) = this(i, j)
        augmented(i, j + n) = if (i == j) 1 else 0
      })
    })
    // perform Gauss-Jordan
    (0 until n).foreach(k => {
      val pivot = augmented(k, k)
      (0 until 2 * n).foreach(j => {
        augmented(k, j) = augmented(k, j) / pivot
      })
      (0 until m).foreach(i => {
        if (i != k) {
          val factor = augmented(i, k)
          (0 until 2 * n).foreach(j => {
            augmented(i, j) = augmented(i, j) - factor * augmented(k, j)
          })
        }
      })
    })
    // extract the inverse matrix from the augmented matrix
    val inverse = new Matrix(n, n)
    (0 until n).foreach(i => {
      (0 until n).foreach(j => {
        inverse(i, j) = augmented(i, j + n)
    })
    })
    inverse
  end inverse
end Matrix

// Companion object
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 identity(n: Int) =
    val a = new Matrix(n)
    (0 until n).foreach(i => a(i,i) = 1)
    a
end Matrix

/* 
 * This is 'leastSquares.scala'. 
 * See "Task 2" below for your task.
 *
 */


package matrixInverse

/*
 * Task 2:
 * This task will put your matrix inversion routine into action
 * with least squares estimation. 
 * Here we ask you to implement a method that fits a polynomial 
 *
 *   p(x) = a(0) + a(1)x + a(2)x^2 + a(3)x^3 + ... + a(d)x^d 
 *
 * with degree d to n data points 
 *
 *   (x(0),y(0)),(x(1),y(1)),...,(x(n-1),y(n-1))
 *
 * so that the sum of squares of errors to data is minimized.
 * That is, our goal is to minimize the sum of (p(x(i))-y(i))^2 over 
 * i=0,1,...,n-1. 
 *
 * Hint: 
 * If you have implemented the matrix inversion method, essentially
 * all you need to do is to set up two matrices, X and y, where y
 * is an n-by-1 matrix (a vector) with the y-coordinates y(i)
 * and X is an n-by-(d+1) matrix whose ith row is 
 * x(i)^0,x(i)^1,x(i)^2,...,x(i)^d. 
 * Then the polynomial can be found with a few basic matrix
 * operations. See here:
 *
 * http://en.wikipedia.org/wiki/Polynomial_regression
 * http://en.wikipedia.org/wiki/Ordinary_least_squares
 *
 */

object leastSquares:
  def fitPolynomial(d: Int, x: Array[Double], y: Array[Double]) : Array[Double] =
    require(d > 0 && x.length == y.length)
    
    ???
  end fitPolynomial
end leastSquares