There are four $3D$ vectors ($3\times 1$ matrix) $\vec a$, $\vec b$, $\vec c$ and $\vec u$ and a $3\times 3$ matrix $L$.
Considering the norm of the following cross product: $$ f=\left\lVert ((\vec a-\vec u)-\vec b) \times ((\vec a- \vec u)-L\vec c)\right\rVert_2$$ how can I compute the partial derivatives of $f$ wrt $u$ and $L$, i.e. $\frac{\partial f}{\partial L} $ and $\frac{\partial f}{\partial u} $? One approach I thought of was converting the first component of cross product to a skew symmetric matrix to turn the cross product to a matrix product. But it becomes very complicated to compute. Do you have any idea? Thanks.
Substitute a new variable $\{x=(u-a)\}$, multiply each side of the cross-product by $(-1)$, and square the function to obtain $$\eqalign{ f^2 &= \|(x+b)\times(x+Lc)\|_2^2 \cr &= \|(x\times Lc)+(b\times Lc)+(b\times x)\|_2^2 \cr &= \|XLc+BLc+Bx)\|_2^2 \cr &= (XLc+BLc+Bx):(XLc+BLc+Bx) \cr &= w:w \cr }$$ where I've used $(B,X)$ to represent the skew-symmetric matrices corresponding to cross-products with respect to the vectors $(b,x)$, and a colon to represent the trace product, e.g. $$A:M={\rm tr}(A^TM)$$ Let's find the differential and gradient wrt $L$ $$\eqalign{ 2f\,df &= 2w:dw \cr &= 2w:(X\,dL\,c+B\,dL\,c) \cr &= 2(X^T+B^T)wc^T:dL \cr &= 2(X^T+B^T)(XLc+BLc+Bx)c^T:dL \cr\cr \frac{\partial f}{\partial L} &= (X^T+B^T)(XLc+BLc+Bx)c^T\,f^{-1} \cr &= -(x+b)\times\Big((x\times Lc)+(b\times Lc)+(b\times x)\Big)c^T\,f^{-1} \cr &= \Bigg[\Big(a-u-b\Big)\times\Big(\big((u-a)\times Lc\big)+(b\times Lc)+\big(b\times(u-a)\big)\Big)\Bigg]c^T\,f^{-1} \cr }$$ I'll leave $\frac{\partial f}{\partial u}$ as an exercise for you; it's pretty much the same technique.
Here are a few hints $$\eqalign{ X &= U + A \cr \frac{\partial f}{\partial x} &= \frac{\partial f}{\partial u} \cr }$$
Update (responding to your comment)
The third-order Levi-Civita tensor $(E)$ is very useful for working with cross products.
Here are a few of its properties $$\eqalign{ B &= -E\cdot b \cr Bx &= b\times x = E:bx^T = -b\cdot E\cdot x \cr b &= -\frac{1}{2}B:E \implies B:E = -2b \cr }$$ Since any matrix can be decomposed into its symmetric and skew parts $$\eqalign{ M &= \frac{1}{2}(M+M^T) + \frac{1}{2}(M-M^T) \cr &= {\rm sym}(M) + {\rm skew}(M) \cr &= Y + K \cr }$$ and contracting with $E$ will annihilate the symmetric part $$\eqalign{M:E &= K:E \cr}$$ So now if we set $$M=wc^TL^T$$ we can write the result from your comment as $$\eqalign{ f\,df &= M:dX + B^Tw:dx \cr &= -M:E\,dx + B^Tw:dx \cr &= (-K:E + B^Tw):dx \cr \frac{\partial f}{\partial x} &= (2k + B^Tw)\,f^{-1} \cr }$$ where $k$ is the vector associated with the skew matrix $K={\rm skew}(M)=-E\cdot k$.
I'll leave it to you to convert this expression back into terms of the original problem variables.