### APDL Pseudoinverse Least Square Fit

Curve fitting is a frequently used tool in engineering. I wished my linear algebra teacher taught me pseudoinverse. As used in the previous blog post script, it computes the least-square curve fit for linear equations. This can come in handy for fitting 2D or even 9D variables.

For a simple 2D straight line, the equation is $y = m \cdot x + c$. Writing it matrix form, we have:

$y = \begin{bmatrix} x & 1 \end{bmatrix} \cdot \begin{bmatrix} m\\c \end{bmatrix}$

With multiple data points in space (e.g. $(x_1,y_1), (x_2,y_2),\dots$ ) the equation looks like this...

$\begin{bmatrix}y_1\\y_2 \\y_3\\ \vdots \end{bmatrix} = \begin{bmatrix} x_1 & 1\\ x_2 & 1 \\ x_3 & 1 \\ \vdots & \vdots \end{bmatrix} \cdot \begin{bmatrix} m\\c \end{bmatrix}$

This matrix looks similar to $y = A \cdot b$ where our

$A = \begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ x_3 & 1 \\ \vdots & \vdots \end{bmatrix}$ and $b = \begin{bmatrix} m\\c \end{bmatrix}$

If $A$ is a square matrix, we can inverse $A$ to get $A^{-1} \cdot y = b$. However when the $A$ matrix is not square, it is over constrained. There are more equations than unknowns like the problem at hand above. That's where the Moore-Penrose inverse does it's magic. From Wikipedia, the pseudoinverse is written as:

$A^{+} = (A^{T} \cdot A)^{-1} \cdot A^{T}$

where $A^{+} \cdot A = I$

Rewriting the linear equation to solve for $b$:

$(A^{T} \cdot A)^{-1} \cdot A^{T} \cdot y = b$

To program it in Ansys efficiently, the equation is rearranged as:

$(A^{T} \cdot A) \cdot b\ = A^{T} \cdot y$

APDL script
Here's an APDL script that works through it.

!! Known answers to shoot for
!! y = mm*x+cc
mm = 12
cc = 34

! Setup A Matrix and fake data with known answers
/prep7
nval = 20
*dim, Amat,, nval, 2
*dim, y,, nval, 1
*do,ct,1,nval
Amat(ct,1) = ct ! x variable
Amat(ct,2) = 1
y(ct,1) = (mm*ct+cc)*(1+0.01*rand(-1,1)) ! adds random noise
*enddo

! Calculates Transpose
*dim, AmatT,, 2, nval
*mfun,AmatT(1,1),tran, Amat(1,1)

! Computes LHS Matrix
*dim, lhs,, 2,2
*moper, lhs(1,1), AmatT(1,1), mult, Amat(1,1) ! LHS = $(A^{T} \cdot A)$

! Computes RHS Matrix
*dim, rhs,, nval,1
*moper, rhs(1,1), AmatT(1,1), mult, y(1,1) ! RHS = $A^{T} \cdot y$

! Solves matrix
*dim, coef,, 2,1
*moper, coef(1,1), lhs(1,1), solv, rhs(1,1) ! LHS*coef = RHS, solve for coef

! Prints out Gradient & Intercept
*stat, coef

Resources
Here's the text file of the above script: [Link]
1. 