I'm trying to implement an algorithm for Convolutional Blind Source Separation (CBSS) based on the ALS algorithm for common BSS on this paper.
On this paper, the problem is formulated by (noise suppressed):
$$\boldsymbol Y=\boldsymbol{AS}$$
which $\boldsymbol A$ and $\boldsymbol S$ must be estimated. The function to be optimized is:
$$J(\boldsymbol Y||\boldsymbol{AS})=\frac{1}{2}||\boldsymbol Y-\boldsymbol{AS}||^2_F+\alpha_A||\boldsymbol A||_{L1}+\alpha_S||\boldsymbol S||_{L1}$$
and the gradients are calculated as:
$$\frac{\partial J(\boldsymbol Y||\boldsymbol{AS})}{\partial \boldsymbol A}=(\boldsymbol{AS}-\boldsymbol Y)\boldsymbol S^T+\alpha_A sign(\boldsymbol A)$$
$$\frac{\partial J(\boldsymbol Y||\boldsymbol{AS})}{\partial \boldsymbol S}=\boldsymbol A^T(\boldsymbol{AS}-\boldsymbol Y)+\alpha_S sign(\boldsymbol S)$$
$S$ is the sources matrix ($N\times P$), $Y$ the observed signals matrix ($M\times P$), $A$ is the mixture matrix ($M\times N$) and the $sign(.)$ function is applied element-wise. $\alpha_A$ and $\alpha_S$ are lagrange multipliers for sparcity regularization on $\boldsymbol A$ and $\boldsymbol S$, respectively.
I want to expand this method for CBSS with $N$ sources and $M$ observed signals and $P$ samples in each signal, where the problem is formulated by ($\circledast$ for circular convolution):
$$\boldsymbol Y=\boldsymbol H\circledast\boldsymbol S$$
where $\boldsymbol H$ is a matrix of discrete LTI systems:
$$\boldsymbol H= \begin{bmatrix} h_{00}[n]&h_{01}[n]&\dots&h_{0M}[n]\\ h_{10}[n]&h_{11}[n]&\dots&h_{1M}[n]\\ \vdots&\vdots&\ddots&\vdots\\ h_{N0}[n]&h_{N1}[n]&\dots&h_{NM}[n]\\ \end{bmatrix}$$
So, the $i$-th row of $\boldsymbol Y$ is calculated by:
$$\boldsymbol y_i=\sum_{j=0}^{N-1}h_{ij}[n]\circledast\boldsymbol s_j[n]=\sum_{j=0}^{N-1}\sum_{p=0}^{P-1}h_{ij}[p]\boldsymbol s_j[n-p]$$
$\boldsymbol s_j[n]$ been the source signal represented by the $j$-th row of $\boldsymbol S$.
So, what should be the gradient functions to expand the ALS algorithm to CBSS in this case?
$$\frac{\partial J(\boldsymbol Y||\boldsymbol H\circledast\boldsymbol S)}{\partial \boldsymbol H}=?$$
$$\frac{\partial J(\boldsymbol Y||\boldsymbol H\circledast\boldsymbol S)}{\partial \boldsymbol S}=?$$
$h_{ij}[n]$ can be seen as a vector of $P$ elements, so $\boldsymbol H$ can be a $M\times N\times P$ tensor
If needed, gradient in $\boldsymbol H$ can be divided in several derivatives of each LTI system, maybe it makes algorithm slower, but feasible:
$$\frac{\partial J(\boldsymbol Y||\boldsymbol H\circledast\boldsymbol S)}{\partial h_{ij}}=?$$
Or even in terms of the gradient of each coefficient of each LTI system:
$$\frac{\partial J(\boldsymbol Y||\boldsymbol H\circledast\boldsymbol S)}{\partial h_{ij}[k]}=?$$
I find a way to estimate the system and the source signals in a convolutional BSS model by using it's Fourier Transform so the convolution becomes a product and the model is estimated in the frequency domain. But is there a way to to it directly from the time domain?
