Skip to main content

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]
Short lecture on it on YouTube: [Link]

Comments

Post a Comment