Deconvolution Question on Article "Deriving Intrinsic Images from Image Sequences" by Yair Weiss there are n derivative filters: $f_i$, and denote $f_i^r$ as $f_i$'s reverse filter such that

$$f_i(x,y)=f_i^r(-x, -y)$$

$r_i, f_i$ given, to find $r$ from the equations:
$$f_i * r = r_i, (1 \leq i \leq n)$$

Professor Weiss recovers $r$ using the following method:
$$r = g * \left(\sum_{i=1}^n f_i^r*r_i\right)$$ and $g$ satisfies:
$$\delta = g * \left(\sum_{i=1}^n f_i^r*f_i\right)$$


Why can the method solve the equation? What books should I refer to?
I want to know how to find the $g$?


Professor Weiss uses these two derivative filters: $[0\ 1 -1]$ and $[0; 1; -1]$.
and in his code the function invDel() returns $g$, which I could not understand. Article link, implementation link.

Professor Weiss's MATLAB code for getting $g$:

function [invK]=invDel2(isize)    % isize is 2 * max(imgWidth, imgHeight)
    K=zeros(isize);
    K(isize/2,isize/2)=-4;
    K(isize/2+1,isize/2)=1;
    K(isize/2,isize/2+1)=1;
    K(isize/2-1,isize/2)=1;
    K(isize/2,isize/2-1)=1;

    Khat=fft2(K);
    I=find(Khat==0);
    Khat(I)=1;
    invKhat=1./Khat;
    invKhat(I)=0;
    invK=ifft2(invKhat);
    invK=-real(invK);
    invK=conv2(invK,[1 0 0;0 0 0;0 0 0],'same');% shift by one
I didn't have time to fully read the article yet I'll try to answer according to your question.

If you take the first equation and convolove it on the left by $ {f}_{i} $ and use the second identity you can see this specific $ g $ holds the equation.

Regarding how to get this $ g $, I think he uses "Deconvolution" in the Fourier domain.