natas_hw
/
openvino
0
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

[Spec][Ref][CPU] Improve Inverse testing suite, add LU pseudocode to specification (#22834) 

### Details:
 - Improved Inverse testing suite for dynamic shapes in the CPU plugin
 - Introduced random tests to PT FE Inerse testing suite
 - Added LU decomposition pseudocode to spec
 - Removed numpy import from opset/opset14.py
- Moved reference inverse outside inverse namespace and removed inverse
namespace
 - Added label propagation tests
 - Changed bf16 threshold restriction
 - Removed dynamic rank check as it's already present in CPU plugin
 - Simplified check for B => P[i] == j
- Did NOT use std::copy_n as ov::bfloat16 is missing difference_value to
use in std implementations

### Tickets:
 - 133090

---------

Co-authored-by: Roman Kazantsev <roman.kazantsev@intel.com>
Co-authored-by: Pawel Raasz <pawel.raasz@intel.com>
Co-authored-by: Michal Lukaszewski <michal.lukaszewski@intel.com>
This commit is contained in:
Piotr Krzemiński 2024-03-28 14:07:57 +01:00 committed by GitHub
parent fba93854cb
commit cb7212b834
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
14 changed files with 499 additions and 391 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
docs/articles_en/documentation/openvino-ir-format/operation-sets/available-opsets/opset14.rst
View File

@ -103,7 +103,7 @@ Table of Contents
* :doc:`I420toRGB <../operation-specs/image/i420-to-rgb-8>`
* :doc:`If <../operation-specs/condition/if-8>`
* :doc:`Interpolate <../operation-specs/image/interpolate-11>`
* :doc:`Inverse <../operation-specs/matrix/Inverse_14>`
* :doc:`Inverse <../operation-specs/matrix/inverse-14>`
* :doc:`IRDFT <../operation-specs/signals/irdft-9>`
* :doc:`IsInf <../operation-specs/comparison/isinf-10>`
* :doc:`IsNaN <../operation-specs/comparison/isnan-10>`

135
docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/Inverse_14.rst
View File

@ -1,135 +0,0 @@
.. {#openvino_docs_ops_matrix_Inverse_14}
Inverse
=======
.. meta::
:description: Learn about Inverse-14 - a matrix operation that computes the inverse or adjoint for one matrix or a batch of input matrice.
**Versioned name**: *Inverse-14*
**Category**: *Matrix*
**Short description**: *Inverse* operation computes either the inverse or adjoint (conjugate transposes) of one or a batch of square invertible matrices.
**Detailed description**: *Inverse* operation computes the inverse of a square matrix. The operation uses LU decomposition with partial pivoting to compute the inverses. (`LU decomposotion with partial pivoting algorithm <https://arxiv.org/abs/2304.03068>`__)
The inverse matrix A^(-1) of a square matrix A is defined as:
.. math::
A \cdot A^{-1} = A^{-1} \cdot A = I
where **I** is the n-dimensional identity matrix.
The inverse matrix exists if and only if the input matrix is invertible. In that case, the inverse is unique. If the matrix is not invertible however, the operation may raise an exception or return an undefined result.
This operation can be used to compute the adjugate matrix instead of the inverse.
The adjugate matrix adj(A) of a square matrix A is defined as:
.. math::
adj(A) = det(A) \cdot A^{-1}
where **A^{-1}** is the matrix inverse of A, and **det(A)** is the determinant of A.
The adjugate matrix exists if and only if the inverse matrix exists.
**Attribute**:
* *adjoint*
* **Description**: Modifies the return value of the operation. If true, the operation returns the adjoint (conjugate transpose) of the input matrices instead of finding the inverse.
* **Range of values**: `true`, `false`
* ``true`` - output adjugate matrix.
* ``false`` - output inverse matrix.
* **Type**: `bool`
* **Default value**: `false`
* **Required**: *No*
**Input**:
* **1**: `input` - A tensor of shape [B1, B2, ..., Bn, ROW, COL] and type `T` representing the input square matrices. **The inner-most 2 dimensions form square matrices and must be of the same size.** B1, B2, ..., Bn represent any amount of batch dimensions (can be 0 for a single matrix input). **Required.**
**Output**:
* **1**: `output` - A tensor with the same type `T` as the input and same shape [B1, B2, ..., Bn, ROW, COL] as the input, representing the inverse matrices (or adjugate matrices) of the input matrices.
**Types**
* **T**: any supported floating-point type.
*Example 1: 2D input matrix.*
.. code-block:: xml
:force:
<layer ... name="Inverse" type="Inverse">
<data/>
<input>
<port id="0" precision="FP32">
<dim>3</dim> <!-- 3 rows of square matrix -->
<dim>3</dim> <!-- 3 columns of square matrix -->
</port>
</input>
<output>
<port id="1" precision="FP32" names="Inverse:0">
<dim>3</dim> <!-- 3 rows of square matrix -->
<dim>3</dim> <!-- 3 columns of square matrix -->
</port>
</output>
</layer>
*Example 2: 3D input tensor with one batch dimension and adjoint=true.*
.. code-block:: xml
:force:
<layer ... name="Inverse" type="Inverse">
<data adjoint="true"/>
<input>
<port id="0" precision="FP32">
<dim>2</dim> <!-- batch size of 2 -->
<dim>4</dim> <!-- 4 rows of square matrix -->
<dim>4</dim> <!-- 4 columns of square matrix -->
</port>
</input>
<output>
<port id="1" precision="FP32" names="Inverse:0">
<dim>2</dim> <!-- batch size of 2 -->
<dim>4</dim> <!-- 4 rows of square matrix -->
<dim>4</dim> <!-- 4 columns of square matrix -->
</port>
</output>
</layer>
*Example 3: 5D input tensor with three batch dimensions.*
.. code-block:: xml
:force:
<layer ... name="Inverse" type="Inverse">
<data/>
<input>
<port id="0" precision="FP32">
<dim>5</dim> <!-- batch size of 5 -->
<dim>4</dim> <!-- batch size of 4 -->
<dim>3</dim> <!-- batch size of 3 -->
<dim>2</dim> <!-- 2 rows of square matrix -->
<dim>2</dim> <!-- 2 columns of square matrix -->
</port>
</input>
<output>
<port id="1" precision="FP32" names="Inverse:0">
<dim>5</dim> <!-- batch size of 5 -->
<dim>4</dim> <!-- batch size of 4 -->
<dim>3</dim> <!-- batch size of 3 -->
<dim>2</dim> <!-- 2 rows of square matrix -->
<dim>2</dim> <!-- 2 columns of square matrix -->
</port>
</output>
</layer>

200
docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/inverse-14.rst Normal file
View File

@ -0,0 +1,200 @@
.. {#openvino_docs_ops_matrix_Inverse_14}
Inverse
=======
.. meta::
:description: Learn about Inverse-14 - a matrix operation that computes the inverse or adjoint for one matrix or a batch of input matrice.
**Versioned name**: *Inverse-14*
**Category**: *Matrix*
**Short description**: *Inverse* operation computes either the inverse or adjoint (conjugate transposes) of one or a batch of square invertible matrices.
**Detailed description**: *Inverse* operation computes the inverse of a square matrix.
The inverse matrix :math:`A^{-1}` of a square n-dimensional matrix :math:`A` is defined as:
.. math::
A \cdot A^{-1} = A^{-1} \cdot A = I
where :math:`I` is the n-dimensional identity matrix.
The inverse matrix exists if and only if the input matrix is invertible (satisfies any of the properties of the *Invertible Matrix Theorem*). In that case, the inverse exists, and is unique. However, if the matrix is not invertible, the operation may raise an exception or return an undefined result. The operation may return slightly different results for the same input data on different devices due to parallelism and different data types implementations.
This operation can be used to compute the adjugate matrix instead of the inverse.
The adjugate matrix :math:`adj(A)` of a square matrix :math:`A` is defined as:
.. math::
adj(A) = det(A) \cdot A^{-1}
where :math:`A^{-1}` is the matrix inverse of :math:`A`, and :math:`det(A)` is the determinant of :math:`A`.
The adjugate matrix exists if and only if the inverse matrix exists.
**Algorithm formulation**:
This operation uses LU decomposition with partial pivoting to compute the inverse (or adjugate) matrix.
.. note::
LU decomposition decomposes matrix :math:`A` into 2 matrices :math:`L` and :math:`U`, where :math:`L` is the lower triangular matrix with all diagonal elements equal to 1 and :math:`U` is the upper triangular matrix.
.. math::
A = L \cdot U
.. note::
LU decomposition allows to easily obtain determinant of :math:`A`. Notice that since :math:`L`` is a lower triangular matrix with all diagonal elements equal to 1, :math:`det(L) = 1`.
.. math::
det(A) = det(L) * det(U) = 1 * det(U) = det(U)
.. note::
To compute the inverse of :math:`A`, given its LU decomposition, it is enough to solve a set of n simple linear equations.
A simple linear equation is of the form :math:`Ax=b`, where x and b are vectors, and x is the solution of the equation. The goal is to compute x, such that it is the i-th column of matrix :math:`A^{-1}`. To do so, let b be a vector of zeros, except for i-th spot that has a value of one (Then b is an i-th column of matrix I).
It is easy to notice that the x-vectors create the matrix :math:`A^{-1}`, and the b-vectors create the Identity matrix.
.. math::
A \cdot A^{-1} = I \implies A^{-1} = [x_1 \& x_2 \& ... \& x_n], A \cdot x_i = b_i, b_i \in col_i(I)
.. note::
Using the LU decomposition, the simple linear equation can be replaced by two, even simpler linear equations.
.. math::
A \cdot x = b \iff L \cdot U \cdot x = b \iff L \cdot y = b, U \cdot x = y
Algorithm pseudocode:
1. Start with original matrix :math:`A`. If the data type of :math:`A` is not f32, convert them to f32 to avoid accumulating rounding errors.
2. Copy initial matrix into matrix :math:`U`. Initialize matrix :math:`L` to be the Identity matrix (zero matrix with all diagonal elements set to 1).
3. Perform LU decomposition with partial pivoting.
* Repeat this step for each column in the input matrix.
* Let *c* be the index of the currently processed column.
* Find the index of the row with the highest value in a given column - *pivot*.
* If :math:`pivot \neq c`, swap the *pivot* and *c* row in :math:`L`. Repeat for :math:`U`. Note that this operation flips the sign of the determinant, so this has to be accounted for.
* Perform standard Gaussian elimination.
4. To obtain the inverse, solve for each column of :math:`A^{-1}` as explained above.
* Solve linear equation :math:`Ly = b` for y (forward substitution)
* Solve linear equation :math:`Ux = y` for x (backward substitution)
* Set x as the corresponding column of the output inverse matrix :math:`A^{-1}`
5. If adjoint == true, then it is necessary to multiply :math:`A^{-1}` by its determinant.
* As explained above, it is enough to compute :math:`det(U)`, since :math:`det(U) = det(A)`.
* :math:`det(U)` is just a multiplication of its diagonal elements.
* Account for each row swap in the LU decomposition step - for every row swap, swap the sign of the determinant.
* Multiply all elements of :math:`A^{-1}` by the determinant to obtain the adjugate matrix.
6. Return the computed matrix. Convert it back from f32 to its original element type.
**Attribute**:
* *adjoint*
* **Description**: Modifies the return value of the operation. If true, the operation returns the adjoint (conjugate transpose) of the input matrices instead of finding the inverse.
* **Range of values**: `true`, `false`
* ``true`` - output adjugate matrix.
* ``false`` - output inverse matrix.
* **Type**: `bool`
* **Default value**: `false`
* **Required**: *No*
**Input**:
* **1**: `input` - A tensor of shape [B1, B2, ..., Bn, ROW, COL] and type `T` representing the input square matrices. **The inner-most 2 dimensions form square matrices and must be of the same size.** B1, B2, ..., Bn represent any amount of batch dimensions (can be 0 for a single matrix input). **Required.**
**Output**:
* **1**: `output` - A tensor with the same type `T` as the input and same shape [B1, B2, ..., Bn, ROW, COL] as the input, representing the inverse matrices (or adjugate matrices) of the input matrices.
**Types**
* **T**: any supported floating-point type. Any type other than f32 will be converted to f32 before executing this op, and then converted back to the original input type to avoid accumulating rounding errors.
*Example 1: 2D input matrix.*
.. code-block:: xml
:force:
<layer ... name="Inverse" type="Inverse">
<data/>
<input>
<port id="0" precision="FP32">
<dim>3</dim> <!-- 3 rows of square matrix -->
<dim>3</dim> <!-- 3 columns of square matrix -->
</port>
</input>
<output>
<port id="1" precision="FP32" names="Inverse:0">
<dim>3</dim> <!-- 3 rows of square matrix -->
<dim>3</dim> <!-- 3 columns of square matrix -->
</port>
</output>
</layer>
*Example 2: 3D input tensor with one batch dimension and adjoint=true.*
.. code-block:: xml
:force:
<layer ... name="Inverse" type="Inverse">
<data adjoint="true"/>
<input>
<port id="0" precision="FP32">
<dim>2</dim> <!-- batch size of 2 -->
<dim>4</dim> <!-- 4 rows of square matrix -->
<dim>4</dim> <!-- 4 columns of square matrix -->
</port>
</input>
<output>
<port id="1" precision="FP32" names="Inverse:0">
<dim>2</dim> <!-- batch size of 2 -->
<dim>4</dim> <!-- 4 rows of square matrix -->
<dim>4</dim> <!-- 4 columns of square matrix -->
</port>
</output>
</layer>
*Example 3: 5D input tensor with three batch dimensions.*
.. code-block:: xml
:force:
<layer ... name="Inverse" type="Inverse">
<data/>
<input>
<port id="0" precision="FP32">
<dim>5</dim> <!-- batch size of 5 -->
<dim>4</dim> <!-- batch size of 4 -->
<dim>3</dim> <!-- batch size of 3 -->
<dim>2</dim> <!-- 2 rows of square matrix -->
<dim>2</dim> <!-- 2 columns of square matrix -->
</port>
</input>
<output>
<port id="1" precision="FP32" names="Inverse:0">
<dim>5</dim> <!-- batch size of 5 -->
<dim>4</dim> <!-- batch size of 4 -->
<dim>3</dim> <!-- batch size of 3 -->
<dim>2</dim> <!-- 2 rows of square matrix -->
<dim>2</dim> <!-- 2 columns of square matrix -->
</port>
</output>
</layer>

134
src/core/reference/include/openvino/reference/inverse.hpp
View File

@ -9,107 +9,30 @@
#include <vector>
#include "openvino/core/shape.hpp"
#include "openvino/reference/convert.hpp"
namespace ov {
namespace reference {
namespace inverse {
namespace internal {
template <typename T>
void lu_decomposition(const T* input,
std::vector<std::vector<T>>& L,
std::vector<std::vector<T>>& U,
std::vector<T>& P,
void lu_decomposition(std::vector<float>& input,
std::vector<float>& L,
std::vector<float>& U,
std::vector<size_t>& P,
bool& sign,
size_t b,
size_t n) {
// Make L identity, U a copy of input and P a range(0, n)
const auto batch_idx = b * n * n;
for (size_t i = 0; i < n; ++i) {
P[i] = static_cast<T>(i);
L[i][i] = T{1};
size_t n,
size_t n_squared);
auto i_idx = i * n;
for (size_t j = 0; j < n; ++j) {
U[i][j] = input[batch_idx + i_idx + j];
}
}
for (size_t k = 0; k < n; ++k) {
// Partial Pivoting
auto pivot_row = k;
for (auto i = k + 1; i < n; ++i) {
if (std::abs(U[i][k]) > std::abs(U[pivot_row][k])) {
pivot_row = i;
}
}
if (pivot_row != k) {
// Swap rows in L, U (A) and P
std::swap(U[k], U[pivot_row]);
std::swap(L[k], L[pivot_row]);
std::swap(P[k], P[pivot_row]);
sign = !sign;
}
for (auto i = k + 1; i < n; ++i) {
L[i][k] = U[i][k] / U[k][k];
for (auto j = k; j < n; ++j) {
U[i][j] -= L[i][k] * U[k][j];
}
}
}
}
template <typename T>
void lu_solve(T* output,
std::vector<std::vector<T>>& L,
std::vector<std::vector<T>>& U,
std::vector<T>& P,
void lu_solve(std::vector<float>& output,
std::vector<float>& L,
std::vector<float>& U,
std::vector<size_t>& P,
size_t b,
size_t n,
size_t column) {
std::vector<T> B(n, T{0});
std::vector<T> X(n, T{0});
std::vector<T> Y(n, T{0});
B[column] = T{1};
size_t n_squared);
// Forward substitution: Ly = Pb
for (size_t i = 0; i < n; ++i) {
Y[i] = B[P[i]];
for (size_t j = 0; j < i; ++j) {
Y[i] -= L[i][j] * Y[j];
}
}
// Backward substitution: Ux = y
for (int i = static_cast<int>(n - 1); i >= 0; --i) {
X[i] = Y[i];
for (size_t j = static_cast<size_t>(i) + 1; j < n; ++j) {
X[i] -= U[i][j] * X[j];
}
X[i] /= U[i][i];
}
size_t batch_idx = b * n * n;
for (size_t row = 0; row < n; ++row) {
output[batch_idx + row * n + column] = X[row];
}
}
template <typename T>
void to_adjoint(T* output, std::vector<std::vector<T>>& U, bool sign, size_t b, size_t n) {
T determinant = sign ? T{1} : T{-1};
for (size_t i = 0; i < n; ++i) {
determinant *= U[i][i];
}
const auto batch_idx = b * n * n;
for (auto idx = batch_idx; idx < batch_idx + n * n; ++idx) {
output[idx] *= determinant;
}
}
void to_adjoint(std::vector<float>& output, std::vector<float>& U, bool sign, size_t b, size_t n, size_t n_squared);
} // namespace internal
/**
@ -124,26 +47,35 @@ void to_adjoint(T* output, std::vector<std::vector<T>>& U, bool sign, size_t b,
template <typename T>
void inverse(const T* input, T* output, const Shape& shape, const bool adjoint) {
const size_t n = shape.back();
const auto total_elements = shape_size<Shape>(shape);
const auto batch_size = total_elements / n / n;
const auto n_squared = n * n;
size_t batch_size = 1;
for (size_t i = 0; i < shape.size() - 2; ++i) {
batch_size = batch_size * shape[i];
}
auto input_conv = std::vector<float>(batch_size * n_squared);
auto output_conv = std::vector<float>(batch_size * n_squared);
convert<T, float>(input, input_conv.data(), batch_size * n_squared);
std::vector<float> L(n_squared);
std::vector<float> U(n_squared);
std::vector<size_t> P(n);
for (size_t b = 0; b < batch_size; ++b) {
std::vector<std::vector<T>> L(n, std::vector<T>(n, T{0}));
std::vector<std::vector<T>> U(n, std::vector<T>(n, T{0}));
std::vector<T> P(n);
bool sign = true;
internal::lu_decomposition(input, L, U, P, sign, b, n);
internal::lu_decomposition(input_conv, L, U, P, sign, b, n, n_squared);
for (size_t column = 0; column < n; ++column) {
internal::lu_solve(output, L, U, P, b, n, column);
}
internal::lu_solve(output_conv, L, U, P, b, n, n_squared);
if (adjoint) {
internal::to_adjoint(output, U, sign, b, n);
internal::to_adjoint(output_conv, U, sign, b, n, n_squared);
}
}
convert<float, T>(output_conv.data(), output, batch_size * n_squared);
}
} // namespace inverse
} // namespace reference
} // namespace ov

124
src/core/reference/src/op/inverse.cpp Normal file
View File

@ -0,0 +1,124 @@
// Copyright (C) 2018-2024 Intel Corporation
// SPDX-License-Identifier: Apache-2.0
//
#include "openvino/reference/inverse.hpp"
namespace ov {
namespace reference {
namespace internal {
void lu_decomposition(std::vector<float>& input,
std::vector<float>& L,
std::vector<float>& U,
std::vector<size_t>& P,
bool& sign,
size_t b,
size_t n,
size_t n_squared) {
// Make L identity, U a copy of input and P a range(0, n)
const auto batch_idx = b * n_squared;
std::fill(L.begin(), L.end(), 0.0f);
memcpy(&U[0], &input[batch_idx], sizeof(float) * n_squared);
for (size_t i = 0; i < n; ++i) {
L[i * n + i] = 1.0f;
P[i] = i;
}
for (size_t k = 0; k < n; ++k) {
// Partial Pivoting
auto pivot_row = k;
auto pivot_idx = pivot_row * n;
const auto k_idx = k * n;
for (auto i = (k + 1) * n, j = k + 1; i < n_squared; i += n, ++j) {
if (std::abs(U[i + k]) > std::abs(U[pivot_idx + k])) {
pivot_row = j;
pivot_idx = pivot_row * n;
}
}
if (pivot_row != k) {
// Swap rows in L, U (A) and P
sign = !sign;
std::swap(P[k], P[pivot_row]);
std::swap_ranges(&U[k_idx], &U[k_idx + n], &U[pivot_idx]);
std::swap_ranges(&L[k_idx], &L[k_idx + n], &L[pivot_idx]);
}
const auto remaining_columns = n - k;
const auto remaining_rows = remaining_columns - 1;
for (size_t i = 0; i < remaining_rows; ++i) {
const auto i_idx = (i + k + 1) * n;
L[i_idx + k] = U[i_idx + k] / U[k_idx + k];
}
for (size_t i = 0; i < remaining_rows * remaining_columns; ++i) {
const auto i_idx = (i / remaining_columns + k + 1) * n;
const auto j_idx = i % remaining_columns + k;
U[i_idx + j_idx] = U[i_idx + j_idx] - L[i_idx + k] * U[k_idx + j_idx];
}
}
}
void lu_solve(std::vector<float>& output,
std::vector<float>& L,
std::vector<float>& U,
std::vector<size_t>& P,
size_t b,
size_t n,
size_t n_squared) {
std::vector<float> X(n);
std::vector<float> Y(n);
for (size_t column = 0; column < n; ++column) {
std::fill(X.begin(), X.end(), 0.0f);
std::fill(Y.begin(), Y.end(), 0.0f);
// Forward substitution: Ly = Pb
for (size_t i = 0; i < n; ++i) {
if (P[i] == column) {
Y[i] = 1.0f;
}
const auto i_idx = i * n;
for (size_t j = 0; j < i; ++j) {
Y[i] -= L[i_idx + j] * Y[j];
}
}
// Backward substitution: Ux = y
for (size_t i = 0; i < n; ++i) {
const auto i_adj = n - i - 1;
const auto i_idx = i_adj * n;
X[i_adj] = Y[i_adj];
for (size_t j = i_adj + 1; j < n; ++j) {
X[i_adj] = X[i_adj] - U[i_idx + j] * X[j];
}
X[i_adj] = X[i_adj] / U[i_idx + i_adj];
}
const auto batch_column_idx = b * n_squared + column;
for (size_t row = 0; row < n; ++row) {
output[batch_column_idx + row * n] = X[row];
}
}
}
void to_adjoint(std::vector<float>& output, std::vector<float>& U, bool sign, size_t b, size_t n, size_t n_squared) {
auto determinant = sign ? 1.0f : -1.0f;
for (size_t i = 0; i < n; ++i) {
determinant *= U[i * n + i];
}
const auto batch_idx = b * n_squared;
for (size_t i = 0; i < n_squared; ++i) {
output[batch_idx + i] *= determinant;
}
}
} // namespace internal
} // namespace reference
} // namespace ov

11
src/core/tests/type_prop/inverse.cpp
View File

@ -33,6 +33,17 @@ TEST_F(TypePropInverseV14Test, input_f64_2x2_constant) {
EXPECT_EQ(op->get_output_partial_shape(0), (ov::PartialShape{2, 2}));
}
TEST_F(TypePropInverseV14Test, symbol_propagation) {
auto input_shape = ov::PartialShape{2, 2};
auto A = std::make_shared<ov::Symbol>(), B = std::make_shared<ov::Symbol>();
set_shape_symbols(input_shape, {A, B});
const auto data = std::make_shared<ov::op::v0::Parameter>(ov::element::f32, input_shape);
const auto op = make_op(data, false);
EXPECT_EQ(op->get_element_type(), ov::element::f32);
EXPECT_EQ(op->get_output_partial_shape(0), (input_shape));
EXPECT_EQ(get_shape_symbols(op->get_output_partial_shape(0)), get_shape_symbols(input_shape));
}
TEST_F(TypePropInverseV14Test, input_f32_static) {
const auto data = std::make_shared<ov::op::v0::Parameter>(ov::element::f32, ov::Shape{10, 4, 4});
const auto op = make_op(data, false);

129
src/plugins/intel_cpu/src/nodes/inverse.cpp
View File

@ -11,9 +11,6 @@
#include "openvino/op/inverse.hpp"
#include "utils/bfloat16.hpp"
// Parallel LU decomposition algorithm with partial pivoting
// Based on the lectures by Prof. Dr. Thomas Huckle, Parallel Numerics
namespace ov {
namespace intel_cpu {
namespace node {
@ -56,7 +53,7 @@ void Inverse::getSupportedDescriptors() {
void Inverse::initSupportedPrimitiveDescriptors() {
m_input_precision = getOriginalInputPrecisionAtPort(INPUT_PORT);
if (!one_of(m_input_precision, ov::element::f32, ov::element::f16, ov::element::bf16)) {
if (m_input_precision != ov::element::f32) {
m_input_precision = ov::element::f32;
}
@ -88,68 +85,56 @@ bool Inverse::created() const {
}
void Inverse::execute(dnnl::stream strm) {
OV_SWITCH(intel_cpu,
InverseExecute,
this,
m_input_precision,
OV_CASE(ov::element::bf16, bfloat16_t),
OV_CASE(ov::element::f16, ov::float16),
OV_CASE(ov::element::f32, float))
inverse();
}
void Inverse::executeDynamicImpl(dnnl::stream strm) {
execute(strm);
}
template <typename T>
void Inverse::inverse() {
const auto* data = getSrcDataAtPortAs<const T>(INPUT_PORT);
auto* output = getDstDataAtPortAs<T>(OUTPUT_PORT);
const auto* data = getSrcDataAtPortAs<const float>(INPUT_PORT);
auto* output = getDstDataAtPortAs<float>(OUTPUT_PORT);
std::vector<T> L(m_side_squared, T{0});
std::vector<T> U(m_side_squared, T{0});
std::vector<T> P(m_side, T{0});
std::vector<float> L(m_side_squared);
std::vector<float> U(m_side_squared);
std::vector<size_t> P(m_side);
for (size_t b = 0; b < m_batches_count; ++b) {
bool sign = true;
lu_decomposition(data, L, U, P, sign, b);
for (size_t column = 0; column < m_side; ++column) {
lu_solve(output, L, U, P, b, column);
}
lu_solve(output, L, U, P, b);
if (m_adjoint) {
// Multiply by det(A) = det(U)
to_adjoint(output, U, sign, b);
}
}
}
template <typename T>
void Inverse::lu_decomposition(const T* data,
std::vector<T>& L,
std::vector<T>& U,
std::vector<T>& P,
void Inverse::lu_decomposition(const float* data,
std::vector<float>& L,
std::vector<float>& U,
std::vector<size_t>& P,
bool& sign,
size_t b) {
// Make L identity, U a copy of data and P a range(0, side)
size_t batch_idx = b * m_side_squared;
const auto batch_idx = b * m_side_squared;
std::fill(L.begin(), L.end(), T{0});
cpu_parallel_memcpy(&U[0], &data[batch_idx], sizeof(T) * m_side_squared);
std::fill(L.begin(), L.end(), 0.0f);
cpu_parallel_memcpy(&U[0], &data[batch_idx], sizeof(float) * m_side_squared);
parallel_for(m_side, [&](size_t i) {
L[i * m_side + i] = T{1};
P[i] = static_cast<T>(i);
L[i * m_side + i] = 1.0f;
P[i] = i;
});
for (size_t k = 0; k < m_side; ++k) {
size_t pivot_row = k;
size_t k_idx = k * m_side;
size_t pivot_idx = pivot_row * m_side;
size_t remaining_columns = m_side - k;
size_t remaining_rows = remaining_columns - 1;
// Partial Pivoting
auto pivot_row = k;
auto pivot_idx = pivot_row * m_side;
const auto k_idx = k * m_side;
// Find maximum value pivot - non-parallel
for (size_t i = (k + 1) * m_side, j = k + 1; i < m_side_squared; i += m_side, ++j) {
@ -169,63 +154,65 @@ void Inverse::lu_decomposition(const T* data,
});
}
const auto remaining_columns = m_side - k;
const auto remaining_rows = remaining_columns - 1;
parallel_for(remaining_rows, [&](size_t i) {
size_t i_idx = (i + k + 1) * m_side;
const auto i_idx = (i + k + 1) * m_side;
L[i_idx + k] = U[i_idx + k] / U[k_idx + k];
});
parallel_for(remaining_rows * remaining_columns, [&](size_t i) {
size_t i_idx = (i / remaining_columns + k + 1) * m_side;
size_t j_idx = i % remaining_columns + k;
const auto i_idx = (i / remaining_columns + k + 1) * m_side;
const auto j_idx = i % remaining_columns + k;
U[i_idx + j_idx] = U[i_idx + j_idx] - L[i_idx + k] * U[k_idx + j_idx];
});
}
}
template <typename T>
void Inverse::lu_solve(T* output, std::vector<T>& L, std::vector<T>& U, std::vector<T>& P, size_t b, size_t column) {
std::vector<T> B(m_side, T{0});
std::vector<T> X(m_side, T{0});
std::vector<T> Y(m_side, T{0});
B[column] = T{1};
void Inverse::lu_solve(float* output, std::vector<float>& L, std::vector<float>& U, std::vector<size_t>& P, size_t b) {
parallel_for(m_side, [&](size_t column) {
std::vector<float> X(m_side, 0.0f);
std::vector<float> Y(m_side, 0.0f);
// Forward substitution: Ly = Pb - not possible to be parallel
for (size_t i = 0; i < m_side; ++i) {
Y[i] = B[P[i]];
size_t i_idx = i * m_side;
for (size_t j = 0; j < i; ++j) {
Y[i] = Y[i] - L[i_idx + j] * Y[j];
// Forward substitution: Ly = Pb
for (size_t i = 0; i < m_side; ++i) {
if (P[i] == column) {
Y[i] = 1.0f;
}
const auto i_idx = i * m_side;
for (size_t j = 0; j < i; ++j) {
Y[i] = Y[i] - L[i_idx + j] * Y[j];
}
}
}
// Backward substitution: Ux = y - not possible to be parallel
for (size_t i = 0; i < m_side; ++i) {
size_t i_adj = m_side - i - 1;
size_t i_idx = i_adj * m_side;
X[i_adj] = Y[i_adj];
for (size_t j = i_adj + 1; j < m_side; ++j) {
X[i_adj] = X[i_adj] - U[i_idx + j] * X[j];
// Backward substitution: Ux = y
for (size_t i = 0; i < m_side; ++i) {
size_t i_adj = m_side - i - 1;
size_t i_idx = i_adj * m_side;
X[i_adj] = Y[i_adj];
for (size_t j = i_adj + 1; j < m_side; ++j) {
X[i_adj] = X[i_adj] - U[i_idx + j] * X[j];
}
X[i_adj] = X[i_adj] / U[i_idx + i_adj];
}
X[i_adj] = X[i_adj] / U[i_idx + i_adj];
}
// Substitute back to get result
size_t batch_idx = b * m_side_squared;
parallel_for(m_side, [&](size_t row) {
output[batch_idx + row * m_side + column] = X[row];
// Substitute back to get result
const auto batch_column_idx = b * m_side_squared + column;
for (size_t row = 0; row < m_side; ++row) {
output[batch_column_idx + row * m_side] = X[row];
}
});
}
template <typename T>
void Inverse::to_adjoint(T* output, std::vector<T>& U, bool sign, size_t b) {
T determinant = sign ? T{1} : T{-1};
void Inverse::to_adjoint(float* output, std::vector<float>& U, bool sign, size_t b) {
auto determinant = sign ? 1.0f : -1.0f;
for (size_t i = 0; i < m_side; ++i) {
determinant = determinant * U[i * m_side + i];
}
size_t batch_idx = b * m_side_squared;
const auto batch_idx = b * m_side_squared;
parallel_for(m_side_squared, [&](size_t i) {
output[batch_idx + i] = output[batch_idx + i] * determinant;
});

22
src/plugins/intel_cpu/src/nodes/inverse.hpp
View File

@ -48,24 +48,18 @@ private:
size_t m_batches_count = 0;
// Helper functions
template <typename T>
void inverse();
template <typename T>
void lu_decomposition(const T* data, std::vector<T>& L, std::vector<T>& U, std::vector<T>& P, bool& sign, size_t b);
void lu_decomposition(const float* data,
std::vector<float>& L,
std::vector<float>& U,
std::vector<size_t>& P,
bool& sign,
size_t b);
template <typename T>
void to_adjoint(T* output, std::vector<T>& U, bool sign, size_t b);
void to_adjoint(float* output, std::vector<float>& U, bool sign, size_t b);
template <typename T>
void lu_solve(T* output, std::vector<T>& L, std::vector<T>& U, std::vector<T>& P, size_t b, size_t column);
template <typename T>
struct InverseExecute {
void operator()(Inverse* node) {
node->inverse<T>();
}
};
void lu_solve(float* output, std::vector<float>& L, std::vector<float>& U, std::vector<size_t>& P, size_t b);
};
} // namespace node

19
src/plugins/intel_cpu/tests/functional/custom/single_layer_tests/instances/common/inverse.cpp
View File

@ -8,22 +8,23 @@ namespace {
using ov::test::InverseLayerTest;
const auto shapes = testing::Values(ov::Shape{10, 10}, ov::Shape{5, 7, 7}, ov::Shape{5, 4, 3, 3}, ov::Shape{100, 2, 2});
const std::vector<std::vector<ov::test::InputShape>> input_shapes = {
{{{5, 4, 4}, {{5, 4, 4}}}},
{{{20, 3, 3}, {{20, 3, 3}}}},
{{{ov::Dimension{1, 70}, -1}, {{3, 3}, {5, 5}, {4, 4}}}},
{{{-1, ov::Dimension{1, 70}, -1}, {{3, 3, 3}, {4, 4, 4}, {5, 5, 5}}}}};
const auto shapes = testing::ValuesIn(input_shapes);
const auto dtypes = testing::Values(ov::element::f32, ov::element::f16, ov::element::bf16);
const auto adjoint = testing::Values(false, true);
const auto test_static = testing::Values(true);
const auto test_dynamic = testing::Values(false);
const auto seed = testing::Values(0u, 1u, 3u);
const auto seed = testing::Values(1, 2, 3);
const auto device_cpu = testing::Values(ov::test::utils::DEVICE_CPU);
const auto params_static = ::testing::Combine(shapes, dtypes, adjoint, test_static, seed, device_cpu);
const auto params_dynamic = ::testing::Combine(shapes, dtypes, adjoint, test_dynamic, seed, device_cpu);
const auto params = ::testing::Combine(shapes, dtypes, adjoint, seed, device_cpu);
INSTANTIATE_TEST_SUITE_P(smoke_InverseStatic, InverseLayerTest, params_static, InverseLayerTest::getTestCaseName);
INSTANTIATE_TEST_SUITE_P(smoke_InverseDynamic, InverseLayerTest, params_dynamic, InverseLayerTest::getTestCaseName);
INSTANTIATE_TEST_SUITE_P(smoke_InverseStatic, InverseLayerTest, params, InverseLayerTest::getTestCaseName);
} // namespace

4
src/plugins/intel_cpu/tests/functional/shared_tests_instances/skip_tests_config.cpp
View File

@ -197,9 +197,7 @@ std::vector<std::string> disabledTestPatterns() {
R"(.*smoke_LPT/InterpolateTransformation.*)",
// Issue: 129931
R"(smoke_LPT/ConvolutionTransformation.CompareWithRefImpl/f32_\[.*,3,16,16\]_CPU_f32_rank=4D_fq_on_data=\{level=256_shape=\[1\]_input_low=\{ 0 \}_input_high=\{ 255 \}_output_low=\{ .*18.7 \}_output_high\{ 18.8 \}_precision=\}_fq_on_weights=\{_255_\[6,1,1,1\]_\{ .*1.52806e.*39, .*0.2, .*0.3, .*0.3, .*0.2, .*0.1 \}_\{ 1.52806e.*39, 0.2, 0.3, 0.3, 0.2, 0.1 \}\})",
// Issue: 132494
R"(.*smoke_Inverse.*bf16.*)",
// Issue: CVS-133173
// Issue: 133173
R"(.*smoke_ScaledAttn_CPU/ScaledAttnLayerCPUTest.CompareWithRefs/netPRC=bf16.*has_scale=0.*)",
R"(.*smoke_LPT_4D/ConvolutionBackpropDataTransformation.CompareWithRefImpl/f32_\[1,8,16,16\]_CPU_f32_\[16,16\]_level=256_shape=\[.*\]_input_low=\{ 0 \}_input_high=\{ 25.5 \}_output_low=\{ 0 \}_output_high\{ 25.5 \}_precision=__255_\[.*\]_\{ -12.7 \}_\{ 12.7 \}_\{\}.*)",
R"(.*smoke_LPT_4D/ConvolutionBackpropDataTransformation.CompareWithRefImpl/f32_\[1,8,16,16\]_CPU_f32_\[16,16\]_level=256_shape=\[1,1,1,1\]_input_low=\{ 0 \}_input_high=\{ 255 \}_output_low=\{ -12.7 \}_output_high\{ 12.8 \}_precision=.*)",

2
src/plugins/template/backend/ops/inverse.cpp
View File

@ -17,7 +17,7 @@ inline bool evaluate(const std::shared_ptr<ov::op::v14::Inverse>& op,
const auto out_shape = ov::op::v14::shape_infer(op.get(), input_shapes).front().to_shape();
outputs[0].set_shape(out_shape);
ov::reference::inverse::inverse<T>(inputs[0].data<const T>(), outputs[0].data<T>(), out_shape, op->get_adjoint());
ov::reference::inverse<T>(inputs[0].data<const T>(), outputs[0].data<T>(), out_shape, op->get_adjoint());
return true;
}

13
src/tests/functional/shared_test_classes/include/shared_test_classes/single_op/inverse.hpp
View File

@ -13,12 +13,11 @@
namespace ov {
namespace test {
using InverseTestParams = typename std::tuple<ov::Shape, // input shape
ov::element::Type, // element type
bool, // adjoint
bool, // test_static
unsigned int, // seed
std::string // device_name
using InverseTestParams = typename std::tuple<std::vector<InputShape>, // input shape
ov::element::Type, // element type
bool, // adjoint
int32_t, // seed
std::string // device_name
>;
class InverseLayerTest : public testing::WithParamInterface<InverseTestParams>, virtual public SubgraphBaseTest {
@ -31,7 +30,7 @@ protected:
void generate_inputs(const std::vector<ov::Shape>& target_shapes) override;
private:
unsigned int m_seed;
int32_t m_seed;
};
} // namespace test
} // namespace ov

66
src/tests/functional/shared_test_classes/src/single_op/inverse.cpp
View File

@ -7,26 +7,26 @@
#include <numeric>
#include <random>
#include "common_test_utils/ov_tensor_utils.hpp"
using namespace ov::test;
namespace ov {
namespace test {
std::string InverseLayerTest::getTestCaseName(const testing::TestParamInfo<InverseTestParams>& obj) {
ov::Shape input_shape;
std::vector<InputShape> input_shape;
ov::element::Type element_type;
bool adjoint;
bool test_static;
unsigned int seed;
int32_t seed;
std::string device_name;
std::tie(input_shape, element_type, adjoint, test_static, seed, device_name) = obj.param;
std::tie(input_shape, element_type, adjoint, seed, device_name) = obj.param;
const char separator = '_';
std::ostringstream result;
result << "IS=" << input_shape.to_string() << separator;
result << "IS=" << input_shape[0].first.to_string() << separator;
result << "dtype=" << element_type.to_string() << separator;
result << "adjoint=" << ov::test::utils::bool2str(adjoint) << separator;
result << "static=" << ov::test::utils::bool2str(test_static) << separator;
result << "seed=" << seed << separator;
result << "device=" << device_name;
@ -34,68 +34,36 @@ std::string InverseLayerTest::getTestCaseName(const testing::TestParamInfo<Inver
}
void InverseLayerTest::SetUp() {
InverseTestParams test_params;
ov::Shape input_shape;
std::vector<InputShape> input_shape;
ov::element::Type element_type;
bool adjoint;
bool test_static;
std::tie(input_shape, element_type, adjoint, test_static, m_seed, targetDevice) = GetParam();
std::tie(input_shape, element_type, adjoint, m_seed, targetDevice) = GetParam();
std::vector<InputShape> in_shapes;
if (!test_static) {
in_shapes.push_back({{}, {input_shape}});
} else {
in_shapes.push_back({input_shape, {input_shape}});
}
init_input_shapes(in_shapes);
init_input_shapes(input_shape);
const auto data = std::make_shared<ov::op::v0::Parameter>(element_type, input_shape);
const auto data = std::make_shared<ov::op::v0::Parameter>(element_type, input_shape[0].first);
data->set_friendly_name("data");
auto inverse = std::make_shared<ov::op::v14::Inverse>(data, adjoint);
const auto inverse = std::make_shared<ov::op::v14::Inverse>(data, adjoint);
ov::ResultVector results{std::make_shared<ov::op::v0::Result>(inverse)};
function = std::make_shared<ov::Model>(results, ParameterVector{data}, "InverseTestCPU");
}
template <typename T>
void fill_data(T* dst, const size_t count, const unsigned int seed) {
std::mt19937 gen(seed);
std::uniform_real_distribution<float> dis(0.0f, 10.0f);
for (size_t i = 0; i < count; i++) {
dst[i] = static_cast<T>(dis(gen));
}
}
void InverseLayerTest::generate_inputs(const std::vector<ov::Shape>& targetInputStaticShapes) {
inputs.clear();
const auto& func_inputs = function->inputs();
const auto& func_input = func_inputs[0];
const auto& name = func_input.get_node()->get_friendly_name();
const auto& in_prc = func_input.get_element_type();
auto tensor = ov::Tensor(in_prc, targetInputStaticShapes[0]);
size_t count = std::accumulate(targetInputStaticShapes[0].begin(),
targetInputStaticShapes[0].end(),
1,
std::multiplies<size_t>());
switch (in_prc) {
case ov::element::f32:
fill_data(tensor.data<ov::element_type_traits<ov::element::f32>::value_type>(), count, m_seed);
break;
case ov::element::f16:
fill_data(tensor.data<ov::element_type_traits<ov::element::f16>::value_type>(), count, m_seed);
break;
case ov::element::bf16:
fill_data(tensor.data<ov::element_type_traits<ov::element::bf16>::value_type>(), count, m_seed);
break;
default:
OPENVINO_THROW("Inverse does not support precision ", in_prc, " for the 'data' input.");
}
ov::test::utils::InputGenerateData in_data;
in_data.start_from = 5;
in_data.range = 5;
in_data.resolution = 1;
in_data.seed = m_seed;
auto tensor = ov::test::utils::create_and_fill_tensor(in_prc, targetInputStaticShapes[0], in_data);
inputs.insert({func_input.get_node_shared_ptr(), tensor});
}
} // namespace test

29
tests/layer_tests/pytorch_tests/test_inverse.py
View File

@ -77,3 +77,32 @@ class TestInverse(PytorchLayerTest):
self._test(aten_inverse_out(), None, "aten::linalg_inv",
ie_device, precision, ir_version, trace_model=True, freeze_model=False, kwargs_to_prepare_input={"out": out})
@pytest.mark.parametrize("shape", [
(10, 2, 2),
(3, 5, 5),
(10, 10),
(7, 6, 5, 4, 4)
])
@pytest.mark.parametrize("dtype", [
np.float64,
np.float32
])
@pytest.mark.parametrize("seed", [
1, 2, 3
])
@pytest.mark.parametrize("out", [
False,
True
])
@pytest.mark.nightly
@pytest.mark.precommit
def test_inverse(self, shape, dtype, seed, out, ie_device, precision, ir_version):
rng = np.random.default_rng(seed)
self.input_tensor = rng.uniform(-10.0, 10.0, shape).astype(dtype)
if not out:
self._test(aten_inverse(), None, "aten::linalg_inv",
ie_device, precision, ir_version, trace_model=True, freeze_model=False)
else:
self._test(aten_inverse_out(), None, "aten::linalg_inv",
ie_device, precision, ir_version, trace_model=True, freeze_model=False, kwargs_to_prepare_input={"out": out})