Task

Given a set of data vectors in the following format:

   
  
    
      
        y
        =
        {
        
          y
          
            1
          
        
        ,
        
          y
          
            2
          
        
        ,
        .
        .
        .
        ,
        
          y
          
            n
          
        
        }
        
      
    
    {\displaystyle y=\{y_{1},y_{2},...,y_{n}\}\,}
  


   
  
    
      
        
          X
          
            i
          
        
        =
        {
        
          x
          
            i
            1
          
        
        ,
        
          x
          
            i
            2
          
        
        ,
        .
        .
        .
        ,
        
          x
          
            i
            n
          
        
        }
        ,
        i
        ∈
        1..
        k
        
      
    
    {\displaystyle X_{i}=\{x_{i1},x_{i2},...,x_{in}\},i\in 1..k\,}
  


Compute the vector 



β
=
{

β

1


,

β

2


,
.
.
.
,

β

k


}


{\displaystyle \beta =\{\beta _{1},\beta _{2},...,\beta _{k}\}}

 using ordinary least squares regression using the following equation:

   
  
    
      
        
          y
          
            j
          
        
        =
        
          Σ
          
            i
          
        
        
          β
          
            i
          
        
        ⋅
        
          x
          
            i
            j
          
        
        ,
        j
        ∈
        1..
        n
      
    
    {\displaystyle y_{j}=\Sigma _{i}\beta _{i}\cdot x_{ij},j\in 1..n}
  


You can assume  y  is given to you as a vector (a one-dimensional array), and  X  is given to you as a two-dimensional array (i.e. matrix).




 
;; Solve a linear system AX=B where A is symmetric and positive definite, so it can be Cholesky decomposed.
(defun linsys (A B)
  (let* ((n (car  (array-dimensions A)))
         (m (cadr (array-dimensions B)))
         (y (make-array n        :element-type 'long-float :initial-element 0.0L0))
         (X (make-array `(,n ,m) :element-type 'long-float :initial-element 0.0L0))
         (L (chol A))) ; A=LL'
 
    (loop for col from 0 to (- m 1) do
       ;; Forward substitution: y = L\B
       (loop for k from 0 to (- n 1)
             do (setf (aref y k)
                      (/ (- (aref B k col)
                            (loop for j from 0 to (- k 1)
                                  sum (* (aref L k j)
                                         (aref y j))))
                         (aref L k k))))
 
       ;; Back substitution. x=L'\y
       (loop for k from (- n 1) downto 0
             do (setf (aref X k col)
                      (/ (- (aref y k)
                            (loop for j from (+ k 1) to (- n 1)
                                  sum (* (aref L j k)
                                         (aref X j col))))
                         (aref L k k)))))
    X))
 
;; Solve a linear least squares problem. Ax=b, with A being mxn, with m>n.
;; Solves the linear system A'Ax=A'b.
(defun lsqr (A b)
  (linsys (mmul (mtp A) A)
          (mmul (mtp A) b)))
 