Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Incorect values for Jacobian when one of the outputs is built with intermediate steps #527

Closed
ioanaif opened this issue Feb 2, 2023 · 1 comment · Fixed by #1121
Closed
Assignees
Milestone

Comments

@ioanaif
Copy link
Contributor

ioanaif commented Feb 2, 2023

Code to reproduce:

int add(int a, int b) { return a + b ;}
void tes1(int a, int b, int* res){ 
    res[0] = a*10 + b*b*9;
    res[0] = add(res[0], 10);
    res[1] = a*a*9 + b*b*10;
}
void tes2(int a, int b, int* res){ 
    res[0] = a*10 + b*b*9;
    res[0] = res[0] + 10;
    res[1] = a*a*9 + b*b*10;
}
auto jacobian = clad::jacobian(test1);
int matr[4];
int res[2] = {0, 0};
jacobian.execute(3, 5, res, matr);

Both versions will output: [64, 190, 0, 0] as the Jacobian result gets summed column wise.

However, if the intermediate step is not present, then the corect result is computed: [10, 90, 54, 100]:

void tes3t(int a, int b, int* res){ 
    res[0] = a*10 + b*b*9 + 10;
    res[1] = a*a*9 + b*b*10;
}

In case of test1 and test2 the jacobianMatrix indexing is at fault:

--------------
{
        int _r6 = 1 * 9;
        int _r7 = _r6 * _t5;
        jacobianMatrix[0UL] += _r7;
        int _r8 = _t6 * _r6;
        jacobianMatrix[0UL] += _r8;
        int _r9 = 1 * 10;
        int _r10 = _r9 * _t7;
        jacobianMatrix[1UL] += _r10;
        int _r11 = _t8 * _r9;
        jacobianMatrix[1UL] += _r11;
    }
    {
        int _jac0 = 0;
        int _jac1 = 0;
        add_pullback(_t4, 10, 1, &_jac0, &_jac1);
        int _r4 = _jac0;
        int _r5 = _jac1;
    }
    {
        int _r0 = 1 * 10;
        jacobianMatrix[0UL] += _r0;
        int _r1 = 1 * 9;
        int _r2 = _r1 * _t2;
        jacobianMatrix[1UL] += _r2;
        int _r3 = _t3 * _r1;
        jacobianMatrix[1UL] += _r3;
    }
--------------

In the case of test3 the Jacobian is generated as:

--------------
{
        int _r4 = 1 * 9;
        int _r5 = _r4 * _t4;
        jacobianMatrix[2UL] += _r5;
        int _r6 = _t5 * _r4;
        jacobianMatrix[2UL] += _r6;
        int _r7 = 1 * 10;
        int _r8 = _r7 * _t6;
        jacobianMatrix[3UL] += _r8;
        int _r9 = _t7 * _r7;
        jacobianMatrix[3UL] += _r9;
    }
    {
        int _r0 = 1 * 10;
        jacobianMatrix[0UL] += _r0;
        int _r1 = 1 * 9;
        int _r2 = _r1 * _t2;
        jacobianMatrix[1UL] += _r2;
        int _r3 = _t3 * _r1;
        jacobianMatrix[1UL] += _r3;
    }
--------------

@davidlange6
Copy link
Collaborator

same problem as #473?

@vgvassilev vgvassilev added this to the v1.6 milestone May 21, 2024
@vgvassilev vgvassilev assigned PetroZarytskyi and unassigned vaithak Jul 11, 2024
@vgvassilev vgvassilev modified the milestones: v1.6, v1.7 Jul 11, 2024
@vgvassilev vgvassilev modified the milestones: v1.7, v1.8 Aug 4, 2024
PetroZarytskyi added a commit to PetroZarytskyi/clad that referenced this issue Nov 14, 2024
PetroZarytskyi added a commit to PetroZarytskyi/clad that referenced this issue Nov 19, 2024
PetroZarytskyi added a commit to PetroZarytskyi/clad that referenced this issue Nov 19, 2024
 Previously, jacobians were based on the non-vectorized reverse mode, which was mostly incapable of capturing multiple outputs. The implementation worked in a few particular cases. For example, it was not possible to differentiate function calls or declare variables inside the original function body.
 This PR implements jacobians using the vectorized forward mode. At the very least, this will solve the issues described above and give a way forward to solve other ones. This also means introducing features to the vectorized fwd mode will introduce the same features to jacobians.
 Let's take a look at the new signature of jacobians. Since the function to be differentiated is expected to have multiple outputs, we should expect the output in the form of array/pointer/reference parameters (just like before). And for every output parameter, we should generate a corresponding adjoint parameter for the user to acquire the results. Since there is no way to specify which parameters are used as output and which are not, adjoints are generated for all array/pointer/reference parameters. For example:

```
void f(double a, double b, double* c)
 -->
void f_jac(double a, double b, double* c, <matrix<double>* _d_c)
```

or

```
void f(double a, double b, double* c, double[7] t)
 -->
void f_jac(double a, double b, double* c, double[7] t,
 array_ref<matrix<double>> _d_c, matrix<double>* _d_t)
```

This signature is also similar to the old one. e.g.
```
df.execute(a, b, c, result); // old behavior
df.execute(a, b, c, &result); // new behavior
```
However, the behavior differs for multiple output parameters, which the old jacobians did not support.

 Note: the same functionality can be achieved by using the vectorized reverse mode, which should probably be implemented at some point. However, the old code for jacobians is unlikely to be useful for that, and there is not much point in keeping it.

Fixes vgvassilev#472, vgvassilev#1057, vgvassilev#480, vgvassilev#527
PetroZarytskyi added a commit to PetroZarytskyi/clad that referenced this issue Nov 19, 2024
PetroZarytskyi added a commit to PetroZarytskyi/clad that referenced this issue Nov 19, 2024
 Previously, jacobians were based on the non-vectorized reverse mode, which was mostly incapable of capturing multiple outputs. The implementation worked in a few particular cases. For example, it was not possible to differentiate function calls or declare variables inside the original function body.
 This PR implements jacobians using the vectorized forward mode. At the very least, this will solve the issues described above and give a way forward to solve other ones. This also means introducing features to the vectorized fwd mode will introduce the same features to jacobians.
 Let's take a look at the new signature of jacobians. Since the function to be differentiated is expected to have multiple outputs, we should expect the output in the form of array/pointer/reference parameters (just like before). And for every output parameter, we should generate a corresponding adjoint parameter for the user to acquire the results. Since there is no way to specify which parameters are used as output and which are not, adjoints are generated for all array/pointer/reference parameters. For example:

```
void f(double a, double b, double* c)
 -->
void f_jac(double a, double b, double* c, <matrix<double>* _d_c)
```

or

```
void f(double a, double b, double* c, double[7] t)
 -->
void f_jac(double a, double b, double* c, double[7] t,
 array_ref<matrix<double>> _d_c, matrix<double>* _d_t)
```

This signature is also similar to the old one. e.g.
```
df.execute(a, b, c, result); // old behavior
df.execute(a, b, c, &result); // new behavior
```
However, the behavior differs for multiple output parameters, which the old jacobians did not support.

 Note: the same functionality can be achieved by using the vectorized reverse mode, which should probably be implemented at some point. However, the old code for jacobians is unlikely to be useful for that, and there is not much point in keeping it.

Fixes vgvassilev#472, Fixes vgvassilev#1057, Fixes vgvassilev#480, Fixes vgvassilev#527
PetroZarytskyi added a commit to PetroZarytskyi/clad that referenced this issue Nov 19, 2024
vgvassilev pushed a commit that referenced this issue Nov 19, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants