diff --git a/include/pineforge/generic_matrix.hpp b/include/pineforge/generic_matrix.hpp
index c4c48c5..f542890 100644
--- a/include/pineforge/generic_matrix.hpp
+++ b/include/pineforge/generic_matrix.hpp
@@ -127,6 +127,11 @@ inline void sort_impl(std::vector<Row>& data, int column, bool ascending) {
 template <typename T>
 class PineGenericMatrix {
     std::vector<std::vector<T>> data_;
+    bool valid_{false};
+
+    void require_valid() const {
+        if (!valid_) throw std::runtime_error("matrix operation on na ID");
+    }
 
 public:
     ~PineGenericMatrix() = default;
@@ -137,6 +142,7 @@ class PineGenericMatrix {
         PineGenericMatrix m;
         m.data_.assign(static_cast<size_t>(rows),
                        std::vector<T>(static_cast<size_t>(cols), init));
+        m.valid_ = true;
         return m;
     }
 
@@ -148,10 +154,12 @@ class PineGenericMatrix {
         PineGenericMatrix m;
         m.data_.assign(static_cast<size_t>(rows),
                        std::vector<T>(static_cast<size_t>(cols), T{}));
+        m.valid_ = true;
         return m;
     }
 
     T get(int row, int col) const {
+        require_valid();
         if (row < 0 || row >= rows())
             throw std::out_of_range("matrix.get: row index out of range");
         if (col < 0 || col >= columns())
@@ -160,6 +168,7 @@ class PineGenericMatrix {
     }
 
     void set(int row, int col, T val) {
+        require_valid();
         if (row < 0 || row >= rows())
             throw std::out_of_range("matrix.set: row index out of range");
         if (col < 0 || col >= columns())
@@ -168,19 +177,28 @@ class PineGenericMatrix {
     }
 
     void fill(T val) {
+        require_valid();
         for (auto& r : data_) std::fill(r.begin(), r.end(), val);
     }
 
-    int rows() const { return static_cast<int>(data_.size()); }
-    int columns() const { return data_.empty() ? 0 : static_cast<int>(data_[0].size()); }
+    int rows() const {
+        require_valid();
+        return static_cast<int>(data_.size());
+    }
+    int columns() const {
+        require_valid();
+        return data_.empty() ? 0 : static_cast<int>(data_[0].size());
+    }
 
     std::vector<T> row(int idx) const {
+        require_valid();
         if (idx < 0 || idx >= rows())
             throw std::out_of_range("matrix.row: row index out of range");
         return data_[static_cast<size_t>(idx)];
     }
 
     std::vector<T> col(int idx) const {
+        require_valid();
         if (idx < 0 || idx >= columns())
             throw std::out_of_range("matrix.col: column index out of range");
         std::vector<T> out;
@@ -192,12 +210,14 @@ class PineGenericMatrix {
     template <typename U = T,
               typename = std::enable_if_t<!std::is_same_v<U, bool>>>
     const std::vector<T>& row_ref(int idx) const {
+        require_valid();
         if (idx < 0 || idx >= rows())
             throw std::out_of_range("matrix.row_ref: row index out of range");
         return data_[static_cast<size_t>(idx)];
     }
 
     void add_row(int idx, const std::vector<T>& values) {
+        require_valid();
         if (idx < 0 || idx > rows())
             throw std::out_of_range("matrix.add_row: row index out of range");
         if (!data_.empty() && values.size() != static_cast<size_t>(columns()))
@@ -207,6 +227,7 @@ class PineGenericMatrix {
     }
 
     void add_col(int idx, const std::vector<T>& values) {
+        require_valid();
         if (data_.empty())
             throw std::logic_error("matrix.add_col on empty matrix: use add_row first");
         if (idx < 0 || idx > columns())
@@ -225,37 +246,44 @@ class PineGenericMatrix {
     }
 
     void remove_row(int idx) {
+        require_valid();
         if (idx < 0 || idx >= rows())
             throw std::out_of_range("matrix.remove_row: row index out of range");
         detail::erase_row(data_, idx);
     }
 
     void remove_col(int idx) {
+        require_valid();
         if (idx < 0 || idx >= columns())
             throw std::out_of_range("matrix.remove_col: column index out of range");
         detail::erase_col(data_, idx);
     }
 
     void swap_rows(int i, int j) {
+        require_valid();
         if (i < 0 || i >= rows() || j < 0 || j >= rows())
             throw std::out_of_range("matrix.swap_rows: row index out of range");
         detail::swap_rows_impl(data_, i, j);
     }
 
     void swap_columns(int i, int j) {
+        require_valid();
         if (i < 0 || i >= columns() || j < 0 || j >= columns())
             throw std::out_of_range("matrix.swap_columns: column index out of range");
         detail::swap_cols_impl(data_, i, j);
     }
 
     [[nodiscard]] PineGenericMatrix copy() const {
+        require_valid();
         PineGenericMatrix m;
         m.data_ = data_;
+        m.valid_ = true;
         return m;
     }
 
     [[nodiscard]] PineGenericMatrix submatrix(int from_row, int to_row,
                                               int from_col, int to_col) const {
+        require_valid();
         if (from_row < 0 || to_row > rows())
             throw std::out_of_range("matrix.submatrix: row index out of range");
         if (from_col < 0 || to_col > columns())
@@ -266,18 +294,23 @@ class PineGenericMatrix {
             throw std::invalid_argument("matrix.submatrix: from_col must be <= to_col");
         PineGenericMatrix m;
         m.data_ = detail::copy_submatrix(data_, from_row, to_row, from_col, to_col);
+        m.valid_ = true;
         return m;
     }
 
     [[nodiscard]] PineGenericMatrix transpose() const {
         static_assert(std::is_default_constructible_v<T>,
                       "matrix.transpose: requires default-constructible element type");
+        require_valid();
         PineGenericMatrix m;
         m.data_ = detail::transpose_impl(data_, rows(), columns(), T{});
+        m.valid_ = true;
         return m;
     }
 
     [[nodiscard]] PineGenericMatrix concat(const PineGenericMatrix& other, bool horizontal) const {
+        require_valid();
+        other.require_valid();
         if (horizontal) {
             if (rows() != other.rows())
                 throw std::invalid_argument("matrix.concat: row count mismatch");
@@ -293,20 +326,27 @@ class PineGenericMatrix {
     void reshape(int new_rows, int new_cols) {
         static_assert(std::is_default_constructible_v<T>,
                       "matrix.reshape: requires default-constructible element type");
+        require_valid();
         detail::reshape_impl(data_, new_rows, new_cols, T{});
     }
 
-    void reverse() { std::reverse(data_.begin(), data_.end()); }
+    void reverse() { require_valid(); std::reverse(data_.begin(), data_.end()); }
 
     void sort(int column, bool ascending = true) {
         static_assert(std::is_same_v<T, int> ||
                       std::is_same_v<T, bool> ||
                       std::is_same_v<T, std::string>,
                       "matrix.sort: requires int, bool, or std::string element type");
+        require_valid();
         detail::sort_impl(data_, column, ascending);
     }
 
-    int elements_count() const { return detail::elements_count_impl(data_); }
+    int elements_count() const {
+        require_valid();
+        return detail::elements_count_impl(data_);
+    }
+
+    [[nodiscard]] bool is_na() const noexcept { return !valid_; }
 };
 
 // PineGenericMatrix<bool> uses std::vector<char> as the row container because
@@ -317,6 +357,11 @@ class PineGenericMatrix {
 template <>
 class PineGenericMatrix<bool> {
     std::vector<std::vector<char>> data_;
+    bool valid_{false};
+
+    void require_valid() const {
+        if (!valid_) throw std::runtime_error("matrix operation on na ID");
+    }
 
 public:
     [[nodiscard]] static PineGenericMatrix new_(int rows, int cols, bool init) {
@@ -325,6 +370,7 @@ class PineGenericMatrix<bool> {
         PineGenericMatrix m;
         m.data_.assign(static_cast<size_t>(rows),
                        std::vector<char>(static_cast<size_t>(cols), init ? 1 : 0));
+        m.valid_ = true;
         return m;
     }
 
@@ -333,6 +379,7 @@ class PineGenericMatrix<bool> {
     }
 
     bool get(int row, int col) const {
+        require_valid();
         if (row < 0 || row >= rows())
             throw std::out_of_range("matrix.get: row index out of range");
         if (col < 0 || col >= columns())
@@ -341,6 +388,7 @@ class PineGenericMatrix<bool> {
     }
 
     void set(int row, int col, bool val) {
+        require_valid();
         if (row < 0 || row >= rows())
             throw std::out_of_range("matrix.set: row index out of range");
         if (col < 0 || col >= columns())
@@ -349,14 +397,22 @@ class PineGenericMatrix<bool> {
     }
 
     void fill(bool val) {
+        require_valid();
         char c = val ? 1 : 0;
         for (auto& r : data_) std::fill(r.begin(), r.end(), c);
     }
 
-    int rows() const { return static_cast<int>(data_.size()); }
-    int columns() const { return data_.empty() ? 0 : static_cast<int>(data_[0].size()); }
+    int rows() const {
+        require_valid();
+        return static_cast<int>(data_.size());
+    }
+    int columns() const {
+        require_valid();
+        return data_.empty() ? 0 : static_cast<int>(data_[0].size());
+    }
 
     std::vector<bool> row(int idx) const {
+        require_valid();
         if (idx < 0 || idx >= rows())
             throw std::out_of_range("matrix.row: row index out of range");
         std::vector<bool> out;
@@ -367,6 +423,7 @@ class PineGenericMatrix<bool> {
     }
 
     std::vector<bool> col(int idx) const {
+        require_valid();
         if (idx < 0 || idx >= columns())
             throw std::out_of_range("matrix.col: column index out of range");
         std::vector<bool> out;
@@ -379,6 +436,7 @@ class PineGenericMatrix<bool> {
     // returning const std::vector<bool>&.
 
     void add_row(int idx, const std::vector<bool>& values) {
+        require_valid();
         if (idx < 0 || idx > rows())
             throw std::out_of_range("matrix.add_row: row index out of range");
         if (!data_.empty() && values.size() != static_cast<size_t>(columns()))
@@ -391,6 +449,7 @@ class PineGenericMatrix<bool> {
     }
 
     void add_col(int idx, const std::vector<bool>& values) {
+        require_valid();
         if (data_.empty())
             throw std::logic_error("matrix.add_col on empty matrix: use add_row first");
         if (idx < 0 || idx > columns())
@@ -408,35 +467,44 @@ class PineGenericMatrix<bool> {
     }
 
     void remove_row(int idx) {
+        require_valid();
         if (idx < 0 || idx >= rows())
             throw std::out_of_range("matrix.remove_row: row index out of range");
         detail::erase_row(data_, idx);
     }
 
     void remove_col(int idx) {
+        require_valid();
         if (idx < 0 || idx >= columns())
             throw std::out_of_range("matrix.remove_col: column index out of range");
         detail::erase_col(data_, idx);
     }
 
     void swap_rows(int i, int j) {
+        require_valid();
         if (i < 0 || i >= rows() || j < 0 || j >= rows())
             throw std::out_of_range("matrix.swap_rows: row index out of range");
         detail::swap_rows_impl(data_, i, j);
     }
 
     void swap_columns(int i, int j) {
+        require_valid();
         if (i < 0 || i >= columns() || j < 0 || j >= columns())
             throw std::out_of_range("matrix.swap_columns: column index out of range");
         detail::swap_cols_impl(data_, i, j);
     }
 
     [[nodiscard]] PineGenericMatrix copy() const {
-        PineGenericMatrix m; m.data_ = data_; return m;
+        require_valid();
+        PineGenericMatrix m;
+        m.data_ = data_;
+        m.valid_ = true;
+        return m;
     }
 
     [[nodiscard]] PineGenericMatrix submatrix(int from_row, int to_row,
                                               int from_col, int to_col) const {
+        require_valid();
         if (from_row < 0 || to_row > rows())
             throw std::out_of_range("matrix.submatrix: row index out of range");
         if (from_col < 0 || to_col > columns())
@@ -447,22 +515,28 @@ class PineGenericMatrix<bool> {
             throw std::invalid_argument("matrix.submatrix: from_col must be <= to_col");
         PineGenericMatrix m;
         m.data_ = detail::copy_submatrix(data_, from_row, to_row, from_col, to_col);
+        m.valid_ = true;
         return m;
     }
 
     void reshape(int new_rows, int new_cols) {
+        require_valid();
         detail::reshape_impl(data_, new_rows, new_cols, static_cast<char>(0));
     }
 
-    void reverse() { std::reverse(data_.begin(), data_.end()); }
+    void reverse() { require_valid(); std::reverse(data_.begin(), data_.end()); }
 
     [[nodiscard]] PineGenericMatrix transpose() const {
+        require_valid();
         PineGenericMatrix m;
         m.data_ = detail::transpose_impl(data_, rows(), columns(), static_cast<char>(0));
+        m.valid_ = true;
         return m;
     }
 
     [[nodiscard]] PineGenericMatrix concat(const PineGenericMatrix& other, bool horizontal) const {
+        require_valid();
+        other.require_valid();
         if (horizontal) {
             if (rows() != other.rows())
                 throw std::invalid_argument("matrix.concat: row count mismatch");
@@ -477,11 +551,22 @@ class PineGenericMatrix<bool> {
 
     template <typename Dummy = void>
     void sort(int /*column*/, bool /*ascending*/ = true) {
+        require_valid();
         static_assert(!std::is_same_v<Dummy, void> && std::is_same_v<Dummy, void>,
                       "matrix.sort: not supported on bool element type");
     }
 
-    int elements_count() const { return detail::elements_count_impl(data_); }
+    int elements_count() const {
+        require_valid();
+        return detail::elements_count_impl(data_);
+    }
+
+    [[nodiscard]] bool is_na() const noexcept { return !valid_; }
 };
 
+template <typename T>
+inline bool is_na(const PineGenericMatrix<T>& matrix) noexcept {
+    return matrix.is_na();
+}
+
 } // namespace pineforge
diff --git a/include/pineforge/matrix.hpp b/include/pineforge/matrix.hpp
index cb91d35..96a2863 100644
--- a/include/pineforge/matrix.hpp
+++ b/include/pineforge/matrix.hpp
@@ -7,6 +7,9 @@ namespace pineforge {
 
 class PineMatrix {
     Eigen::MatrixXd data_;
+    bool valid_{false};
+
+    void require_valid() const;
 
 public:
     // Construction
@@ -67,6 +70,10 @@ class PineMatrix {
     // Count
     int elements_count() const;
 
+    // Pine matrix ID state. A default-constructed value represents ``na``;
+    // matrix.new() returns a valid ID even when its dimensions are 0x0.
+    [[nodiscard]] bool is_na() const noexcept { return !valid_; }
+
     // Properties
     bool is_square() const;
     bool is_identity() const;
@@ -80,8 +87,12 @@ class PineMatrix {
     bool is_zero() const;
 
     // Access to internal data for friend operations
-    const Eigen::MatrixXd& data() const { return data_; }
-    Eigen::MatrixXd& data() { return data_; }
+    const Eigen::MatrixXd& data() const { require_valid(); return data_; }
+    Eigen::MatrixXd& data() { require_valid(); return data_; }
 };
 
+inline bool is_na(const PineMatrix& matrix) noexcept {
+    return matrix.is_na();
+}
+
 } // namespace pineforge
diff --git a/src/matrix.cpp b/src/matrix.cpp
index d4b1032..2351b88 100644
--- a/src/matrix.cpp
+++ b/src/matrix.cpp
@@ -14,37 +14,54 @@ static constexpr double EPS = 1e-10;
 PineMatrix PineMatrix::new_(int rows, int cols, double init_val) {
     PineMatrix m;
     m.data_ = Eigen::MatrixXd::Constant(rows, cols, init_val);
+    m.valid_ = true;
     return m;
 }
 
+void PineMatrix::require_valid() const {
+    if (!valid_) throw std::runtime_error("matrix operation on na ID");
+}
+
 // ── Access ──────────────────────────────────────────────────────────────────
 
 double PineMatrix::get(int row, int col) const {
+    require_valid();
     return data_(row, col);
 }
 
 void PineMatrix::set(int row, int col, double val) {
+    require_valid();
     data_(row, col) = val;
 }
 
 void PineMatrix::fill(double val) {
+    require_valid();
     data_.setConstant(val);
 }
 
 std::vector<double> PineMatrix::row(int idx) const {
+    require_valid();
     Eigen::VectorXd r = data_.row(idx);
     return std::vector<double>(r.data(), r.data() + r.size());
 }
 
 std::vector<double> PineMatrix::col(int idx) const {
+    require_valid();
     Eigen::VectorXd c = data_.col(idx);
     return std::vector<double>(c.data(), c.data() + c.size());
 }
 
 // ── Row/Col ops ─────────────────────────────────────────────────────────────
 
-int PineMatrix::rows() const { return static_cast<int>(data_.rows()); }
-int PineMatrix::columns() const { return static_cast<int>(data_.cols()); }
+int PineMatrix::rows() const {
+    require_valid();
+    return static_cast<int>(data_.rows());
+}
+
+int PineMatrix::columns() const {
+    require_valid();
+    return static_cast<int>(data_.cols());
+}
 
 void PineMatrix::add_row(int idx, const std::vector<double>& values) {
     int r = rows(), c = columns();
@@ -86,24 +103,30 @@ void PineMatrix::remove_col(int idx) {
 // ── Swap ────────────────────────────────────────────────────────────────────
 
 void PineMatrix::swap_rows(int i, int j) {
+    require_valid();
     data_.row(i).swap(data_.row(j));
 }
 
 void PineMatrix::swap_columns(int i, int j) {
+    require_valid();
     data_.col(i).swap(data_.col(j));
 }
 
 // ── Transform ───────────────────────────────────────────────────────────────
 
 PineMatrix PineMatrix::copy() const {
+    require_valid();
     PineMatrix m;
     m.data_ = data_;
+    m.valid_ = true;
     return m;
 }
 
 PineMatrix PineMatrix::submatrix(int from_row, int to_row, int from_col, int to_col) const {
+    require_valid();
     PineMatrix m;
     m.data_ = data_.block(from_row, from_col, to_row - from_row, to_col - from_col);
+    m.valid_ = true;
     return m;
 }
 
@@ -124,13 +147,16 @@ void PineMatrix::reshape(int r, int c) {
 }
 
 void PineMatrix::reverse() {
+    require_valid();
     Eigen::MatrixXd tmp = data_.colwise().reverse();
     data_ = tmp.rowwise().reverse();
 }
 
 PineMatrix PineMatrix::transpose() const {
+    require_valid();
     PineMatrix m;
     m.data_ = data_.transpose();
+    m.valid_ = true;
     return m;
 }
 
@@ -150,6 +176,8 @@ void PineMatrix::sort(int column, bool ascending) {
 }
 
 PineMatrix PineMatrix::concat(const PineMatrix& other, bool horizontal) const {
+    require_valid();
+    other.require_valid();
     PineMatrix m;
     if (horizontal) {
         m.data_.resize(rows(), columns() + other.columns());
@@ -158,17 +186,19 @@ PineMatrix PineMatrix::concat(const PineMatrix& other, bool horizontal) const {
         m.data_.resize(rows() + other.rows(), columns());
         m.data_ << data_, other.data_;
     }
+    m.valid_ = true;
     return m;
 }
 
 // ── Aggregation ─────────────────────────────────────────────────────────────
 
-double PineMatrix::avg() const { return data_.mean(); }
-double PineMatrix::min() const { return data_.minCoeff(); }
-double PineMatrix::max() const { return data_.maxCoeff(); }
-double PineMatrix::sum() const { return data_.sum(); }
+double PineMatrix::avg() const { require_valid(); return data_.mean(); }
+double PineMatrix::min() const { require_valid(); return data_.minCoeff(); }
+double PineMatrix::max() const { require_valid(); return data_.maxCoeff(); }
+double PineMatrix::sum() const { require_valid(); return data_.sum(); }
 
 double PineMatrix::mode() const {
+    require_valid();
     std::map<double, int> freq;
     for (int i = 0; i < rows(); ++i)
         for (int j = 0; j < columns(); ++j)
@@ -195,25 +225,32 @@ bool all_coeffs_finite(const Eigen::MatrixXd& m) {
 // ── Arithmetic ──────────────────────────────────────────────────────────────
 
 PineMatrix PineMatrix::diff(const PineMatrix& other) const {
+    require_valid();
+    other.require_valid();
     PineMatrix m;
     m.data_ = data_ - other.data_;
+    m.valid_ = true;
     return m;
 }
 
 PineMatrix PineMatrix::mult(const PineMatrix& other) const {
+    require_valid();
+    other.require_valid();
     PineMatrix m;
     m.data_ = data_ * other.data_;
+    m.valid_ = true;
     return m;
 }
 
 PineMatrix PineMatrix::pow(int n) const {
     if (!is_square())
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     if (!all_coeffs_finite(data_))
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     // start with identity
     PineMatrix result;
     result.data_ = Eigen::MatrixXd::Identity(rows(), rows());
+    result.valid_ = true;
     for (int i = 0; i < n; ++i)
         result.data_ = result.data_ * data_;
     return result;
@@ -222,6 +259,7 @@ PineMatrix PineMatrix::pow(int n) const {
 // ── Linear algebra ──────────────────────────────────────────────────────────
 
 double PineMatrix::det() const {
+    require_valid();
     if (!is_square())
         return std::numeric_limits<double>::quiet_NaN();
     if (!all_coeffs_finite(data_))
@@ -230,32 +268,38 @@ double PineMatrix::det() const {
 }
 
 PineMatrix PineMatrix::inv() const {
+    require_valid();
     if (!is_square())
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     if (!all_coeffs_finite(data_))
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     PineMatrix m;
     m.data_ = data_.inverse();
+    m.valid_ = true;
     return m;
 }
 
 PineMatrix PineMatrix::pinv() const {
+    require_valid();
     if (!all_coeffs_finite(data_))
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     PineMatrix m;
     m.data_ = data_.completeOrthogonalDecomposition().pseudoInverse();
+    m.valid_ = true;
     return m;
 }
 
 int PineMatrix::rank() const {
+    require_valid();
     if (!all_coeffs_finite(data_))
         return 0;
     return static_cast<int>(Eigen::FullPivLU<Eigen::MatrixXd>(data_).rank());
 }
 
-double PineMatrix::trace() const { return data_.trace(); }
+double PineMatrix::trace() const { require_valid(); return data_.trace(); }
 
 std::vector<double> PineMatrix::eigenvalues() const {
+    require_valid();
     if (!is_square())
         return {};
     if (!all_coeffs_finite(data_))
@@ -285,24 +329,26 @@ std::vector<double> PineMatrix::eigenvalues() const {
 
 PineMatrix PineMatrix::eigenvectors() const {
     if (!is_square())
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     if (!all_coeffs_finite(data_))
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     if (is_symmetric()) {
         Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver(
             data_, Eigen::ComputeEigenvectors);
         if (solver.info() != Eigen::Success)
-            return PineMatrix();
+            return PineMatrix::new_(0, 0, 0.0);
         PineMatrix m;
         m.data_ = solver.eigenvectors();
+        m.valid_ = true;
         return m;
     }
     Eigen::EigenSolver<Eigen::MatrixXd> solver(data_);
     if (solver.info() != Eigen::Success)
-        return PineMatrix();
+        return PineMatrix::new_(0, 0, 0.0);
     auto ev = solver.eigenvectors();
     PineMatrix m;
     m.data_.resize(ev.rows(), ev.cols());
+    m.valid_ = true;
     for (Eigen::Index i = 0; i < ev.rows(); ++i)
         for (Eigen::Index j = 0; j < ev.cols(); ++j)
             m.data_(i, j) = ev(i, j).real();
@@ -312,10 +358,13 @@ PineMatrix PineMatrix::eigenvectors() const {
 // ── Kronecker ───────────────────────────────────────────────────────────────
 
 PineMatrix PineMatrix::kron(const PineMatrix& other) const {
+    require_valid();
+    other.require_valid();
     int ar = rows(), ac = columns();
     int br = other.rows(), bc = other.columns();
     PineMatrix m;
     m.data_.resize(ar * br, ac * bc);
+    m.valid_ = true;
     for (int i = 0; i < ar * br; ++i)
         for (int j = 0; j < ac * bc; ++j)
             m.data_(i, j) = data_(i / br, j / bc) * other.data_(i % br, j % bc);
@@ -325,6 +374,7 @@ PineMatrix PineMatrix::kron(const PineMatrix& other) const {
 // ── Count ───────────────────────────────────────────────────────────────────
 
 int PineMatrix::elements_count() const {
+    require_valid();
     return static_cast<int>(data_.size());
 }
 
@@ -394,6 +444,7 @@ bool PineMatrix::is_binary() const {
 }
 
 bool PineMatrix::is_zero() const {
+    require_valid();
     return data_.isZero(EPS);
 }
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ac7dbeb..86145d4 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -15,6 +15,7 @@ set(TEST_SOURCES
     test_ta_rma_warmup
     test_ta_na_source_rules
     test_matrix
+    test_matrix_na
     test_session_time
     test_time_tradingday
     test_chart_timezone
diff --git a/tests/test_matrix_na.cpp b/tests/test_matrix_na.cpp
new file mode 100644
index 0000000..aef234c
--- /dev/null
+++ b/tests/test_matrix_na.cpp
@@ -0,0 +1,97 @@
+#include <pineforge/generic_matrix.hpp>
+#include <pineforge/matrix.hpp>
+#include <pineforge/na.hpp>
+
+#include <stdexcept>
+#include <string>
+
+using pineforge::PineGenericMatrix;
+using pineforge::PineMatrix;
+using pineforge::is_na;
+
+static void require(bool condition, const char* message) {
+    if (!condition) throw std::runtime_error(message);
+}
+
+template <typename Fn>
+static bool throws_na_id(Fn&& fn) {
+    try {
+        fn();
+    } catch (const std::runtime_error& exc) {
+        return std::string(exc.what()) == "matrix operation on na ID";
+    }
+    return false;
+}
+
+template <typename Matrix, typename Factory>
+static void check_nullable_id_lifecycle(Factory&& make_valid) {
+    Matrix state;
+    require(is_na(state), "default matrix must be a Pine na ID");
+    require(throws_na_id([&] { (void)state.rows(); }),
+            "rows() on na matrix must fail");
+    require(throws_na_id([&] { (void)state.get(-1, 0); }),
+            "na matrix must fail before index validation");
+    require(throws_na_id([&] { state.add_row(-1, {}); }),
+            "na matrix add_row must fail before index validation");
+    require(throws_na_id([&] { state.swap_rows(-1, 0); }),
+            "na matrix swap_rows must fail before index validation");
+    require(throws_na_id([&] { (void)state.submatrix(-1, 0, 0, 0); }),
+            "na matrix submatrix must fail before index validation");
+    require(throws_na_id([&] { (void)state.copy(); }),
+            "copy() on na matrix must fail");
+
+    Matrix copied_null = state;
+    require(is_na(copied_null), "copying a na ID must preserve na");
+
+    state = make_valid(0, 0);
+    require(!is_na(state), "empty matrix must have a valid ID");
+    require(state.rows() == 0, "empty matrix must have zero rows");
+    require(state.columns() == 0, "empty matrix must have zero columns");
+
+    state = make_valid(1, 1);
+    require(!is_na(state), "new matrix must have a valid ID");
+    require(state.rows() == 1, "new matrix must expose its row count");
+    require(state.columns() == 1,
+            "new matrix must expose its column count");
+    require(!is_na(state.copy()), "matrix.copy() must return a valid ID");
+    require(!is_na(state.transpose()),
+            "matrix.transpose() must return a valid ID");
+    require(!is_na(state.concat(make_valid(1, 1), true)),
+            "matrix.concat() must return a valid ID");
+
+    Matrix null_rhs;
+    require(throws_na_id([&] { (void)state.concat(null_rhs, true); }),
+            "matrix.concat() must reject a na RHS ID");
+
+    state = pineforge::na<Matrix>();
+    require(is_na(state), "reassigning na must clear the matrix ID");
+    require(throws_na_id([&] { (void)state.elements_count(); }),
+            "elements_count() on na matrix must fail");
+}
+
+int main() {
+    check_nullable_id_lifecycle<PineMatrix>(
+        [](int rows, int cols) {
+            return PineMatrix::new_(rows, cols, 0.0);
+        });
+    check_nullable_id_lifecycle<PineGenericMatrix<int>>(
+        [](int rows, int cols) {
+            return PineGenericMatrix<int>::new_(rows, cols, 0);
+        });
+    check_nullable_id_lifecycle<PineGenericMatrix<std::string>>(
+        [](int rows, int cols) {
+            return PineGenericMatrix<std::string>::new_(
+                rows, cols, std::string{});
+        });
+    check_nullable_id_lifecycle<PineGenericMatrix<bool>>(
+        [](int rows, int cols) {
+            return PineGenericMatrix<bool>::new_(rows, cols, false);
+        });
+
+    PineMatrix nonsquare = PineMatrix::new_(1, 2, 0.0);
+    PineMatrix empty_inverse = nonsquare.inv();
+    require(!is_na(empty_inverse),
+            "a valid empty failure result must not become a na ID");
+    require(empty_inverse.rows() == 0 && empty_inverse.columns() == 0,
+            "non-square inverse must preserve its valid empty result");
+}