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:

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

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

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:

where \(A^{+} \cdot A = I \)

Rewriting the linear equation to solve for \(b\):

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

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

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))

*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)

! 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

Here's the text file of the above script:

Short lecture on it on YouTube:

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]**
keep it up

ReplyDeleteThank you a lot, you solve my big problem

ReplyDelete