diff --git a/NUMPY_API_SUMMARY.md b/NUMPY_API_SUMMARY.md
new file mode 100644
index 0000000..9de5720
--- /dev/null
+++ b/NUMPY_API_SUMMARY.md
@@ -0,0 +1,221 @@
+# NumPy-Compatible API Layer - Implementation Summary
+
+## Overview
+This document summarizes the implementation of the NumPy-compatible API layer for the Paper framework. The goal was to provide a familiar NumPy interface while leveraging Paper's out-of-core capabilities for datasets larger than memory.
+
+**Current Limitation:** The implementation currently supports **2D matrices only** (NumPy's ndarray supports N dimensions). This is a foundational version focused on the most common use case for matrix operations. Future enhancements may add full N-dimensional support.
+
+## Key Components Implemented
+
+### 1. Core API Module (`paper/numpy_api.py`)
+
+#### ndarray Class
+The main array class that provides NumPy-like interface with lazy evaluation:
+
+**Properties:**
+- `shape` - Tuple representing array dimensions
+- `dtype` - NumPy data type of elements
+- `ndim` - Number of dimensions (always 2 for matrices)
+- `size` - Total number of elements
+- `T` - Transpose property
+
+**Public Methods:**
+- `to_numpy()` - Convert to NumPy array (loads into memory)
+- `compute()` - Execute lazy computation plan
+- `__add__()` - Element-wise addition operator
+- `__mul__()` - Scalar multiplication operator
+- `__matmul__()` - Matrix multiplication operator
+
+**Internal Methods:**
+- `_materialize()` - Internal method for materialization
+- `_from_plan()` - Class method for creating lazy arrays
+
+#### Array Creation Functions
+- `array(data, dtype)` - Create from data (list, tuple, numpy array)
+- `zeros(shape, dtype)` - Create zeros array
+- `ones(shape, dtype)` - Create ones array
+- `eye(n, dtype)` - Create identity matrix
+- `random_rand(shape, dtype)` - Create random array
+
+#### I/O Functions
+- `load(filepath, shape, dtype)` - Load from binary file
+- `save(filepath, array)` - Save to binary file
+
+#### Helper Functions
+- `dot(a, b)` - Matrix multiplication (alias for @)
+- `add(a, b)` - Addition (alias for +)
+- `multiply(a, b)` - Multiplication (alias for *)
+
+### 2. Bug Fix (`paper/plan.py`)
+
+**Issue:** During implementation, we discovered a pre-existing bug in the scalar multiplication logic where matrices from EagerNodes were being closed prematurely, causing segmentation faults on reuse. This bug existed in the original codebase but became apparent when implementing reusable lazy operations in the NumPy API.
+
+**Fix:** Added check to only close intermediate computed results:
+```python
+# Only close TMP if it's not from an EagerNode
+if not isinstance(self.left, EagerNode):
+    TMP.close()
+```
+
+### 3. Test Suite (`tests/test_numpy_api.py`)
+
+Comprehensive testing with **26 tests** covering:
+
+- **Array Creation (9 tests)**
+  - Creation from lists, numpy arrays
+  - Different dtypes support
+  - Special arrays (zeros, ones, eye, random)
+  - Array properties
+  - Public API (to_numpy method)
+
+- **Operations (6 tests)**
+  - Addition
+  - Scalar multiplication (both left and right)
+  - Matrix multiplication
+  - Chained operations
+  - Transpose
+
+- **Error Handling (3 tests)**
+  - Shape mismatch in addition
+  - Dimension mismatch in matmul
+  - Invalid array creation
+
+- **File I/O (3 tests)**
+  - Loading arrays
+  - Saving arrays
+  - Saving lazy arrays
+
+- **Helper Functions (3 tests)**
+  - dot() function
+  - add() function
+  - multiply() function
+
+- **Large Arrays (2 tests)**
+  - Large array creation (1000x1000)
+  - Large matrix multiplication (500x600 @ 600x400)
+
+### 4. Examples
+
+#### Simple Demo (`examples/simple_demo.py`)
+Quick-start example demonstrating basic usage:
+- Array creation
+- Operation chaining
+- Lazy evaluation
+- Result computation
+
+#### Comprehensive Examples (`examples/numpy_api_example.py`)
+9 detailed examples covering:
+1. Basic array creation and operations
+2. Scalar multiplication
+3. Matrix multiplication
+4. Chained operations with lazy evaluation
+5. Array creation functions
+6. Transpose operation
+7. File I/O operations
+8. Large arrays (out-of-core)
+9. NumPy compatibility comparison
+
+### 5. Documentation (`README.md`)
+
+Updated README with:
+- NumPy API overview
+- Quick start guide
+- Key features list
+- Supported operations
+- Examples reference
+
+## Key Features
+
+✅ **Familiar Interface** - NumPy-like syntax for easy migration  
+✅ **Lazy Evaluation** - Build computation plans without immediate execution  
+✅ **Automatic Optimization** - Operator fusion applied automatically  
+✅ **Out-of-Core Support** - Handle datasets larger than memory  
+✅ **Public API** - Clean separation of public and internal methods  
+✅ **Comprehensive Testing** - 62 total tests (26 new + 36 existing)
+
+## Usage Examples
+
+### Basic Usage
+```python
+from paper import numpy_api as pnp
+import numpy as np
+
+# Create arrays
+a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+
+# Build computation plan (lazy)
+c = (a + b) * 2
+
+# Execute and get result
+result = c.compute()
+numpy_result = result.to_numpy()  # Convert to NumPy array
+```
+
+### Matrix Multiplication
+```python
+# Create matrices
+A = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)  # 2x3
+B = pnp.array([[7, 8], [9, 10], [11, 12]], dtype=np.float32)  # 3x2
+
+# Matrix multiplication
+C = A @ B  # Lazy - 2x2 result
+
+# Compute
+result = C.compute()
+print(result.to_numpy())
+```
+
+### Large Arrays
+```python
+# Create large random arrays (out-of-core)
+a = pnp.random_rand((10000, 10000))
+b = pnp.random_rand((10000, 10000))
+
+# Operations don't load entire arrays into memory
+c = (a + b) * 0.5
+
+# Computation uses disk-backed matrices efficiently
+result = c.compute()
+```
+
+## Test Results
+
+All tests passing:
+```
+----------------------------------------------------------------------
+Ran 62 tests in 0.644s
+
+OK
+
+==================================================
+Tests run: 62
+Failures: 0
+Errors: 0
+Skipped: 0
+All tests passed! ✓
+```
+
+## Benefits for Users
+
+1. **Easy Migration**: Minimal code changes required when migrating from NumPy
+2. **Familiar Syntax**: Use the same operations and methods as NumPy
+3. **Scalability**: Handle datasets that don't fit in memory
+4. **Performance**: Automatic optimizations like operator fusion
+5. **Lazy Evaluation**: Build complex computation plans efficiently
+6. **Clean API**: Well-documented public methods with internal implementation hidden
+
+## Future Enhancements
+
+Potential areas for future development:
+- More NumPy operations (subtraction, division, power, etc.)
+- Broadcasting support
+- Slicing operations
+- Reduction operations (sum, mean, max, min)
+- Element-wise functions (sin, cos, exp, log)
+- Lazy transpose implementation
+- Multi-dimensional array support (currently 2D only)
+
+## Conclusion
+
+The NumPy-compatible API layer successfully provides a familiar interface for Paper framework users while maintaining all the benefits of out-of-core matrix operations. The implementation is well-tested, documented, and ready for production use.
diff --git a/README.md b/README.md
index 8bf65cc..e26face 100644
--- a/README.md
+++ b/README.md
@@ -4,6 +4,70 @@ Paper a lightweight Python framework for performing matrix computations on datas
 
 The architecture is inspired by modern data systems and academic research (e.g., PreVision), with a clear separation between the logical plan, the physical execution backend, and an intelligent optimizer.
 
+## NumPy-Compatible API
+
+Paper now includes a **NumPy-compatible API layer** that provides a familiar interface for users migrating from NumPy or other array libraries. This makes it easy to leverage Paper's out-of-core capabilities with minimal code changes.
+
+### Quick Start
+
+```python
+# Import Paper's NumPy-compatible API
+from paper import numpy_api as pnp
+import numpy as np
+
+# Create arrays (similar to NumPy)
+a = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)
+b = pnp.array([[7, 8, 9], [10, 11, 12]], dtype=np.float32)
+
+# Perform operations with lazy evaluation
+c = (a + b) * 2
+
+# Execute the computation plan
+result = c.compute()
+print(result.to_numpy())
+```
+
+### Key Features
+
+- **Familiar NumPy Interface**: Use the same syntax as NumPy for array creation and operations
+- **Lazy Evaluation**: Build computation plans without executing until `.compute()` is called
+- **Automatic Optimization**: Operator fusion and intelligent caching applied automatically
+- **Out-of-Core Support**: Handle datasets larger than memory seamlessly
+- **Matrix Operations**: Support for addition, scalar multiplication, and matrix multiplication (@)
+
+### Supported Operations
+
+**Array Creation:**
+- `pnp.array(data)` - Create array from data
+- `pnp.zeros(shape)` - Create zeros array
+- `pnp.ones(shape)` - Create ones array
+- `pnp.eye(n)` - Create identity matrix
+- `pnp.random_rand(shape)` - Create random array
+
+**Operations:**
+- `a + b` - Element-wise addition
+- `a * scalar` - Scalar multiplication
+- `a @ b` - Matrix multiplication
+- `a.T` - Transpose
+
+**I/O:**
+- `pnp.load(filepath, shape)` - Load array from file
+- `pnp.save(filepath, array)` - Save array to file
+
+### Examples
+
+See `examples/numpy_api_example.py` for comprehensive examples demonstrating:
+- Basic array operations
+- Chained operations with lazy evaluation
+- Matrix multiplication
+- File I/O
+- Large array handling (out-of-core)
+
+Run the examples:
+```bash
+python examples/numpy_api_example.py
+```
+
 ### Architecture
 
 
diff --git a/examples/numpy_api_example.py b/examples/numpy_api_example.py
new file mode 100644
index 0000000..35c5a7b
--- /dev/null
+++ b/examples/numpy_api_example.py
@@ -0,0 +1,272 @@
+#!/usr/bin/env python3
+"""
+NumPy-Compatible API Example
+
+This example demonstrates how to use the Paper framework's NumPy-compatible API
+for out-of-core matrix operations. The API provides a familiar NumPy interface
+while handling datasets larger than memory.
+"""
+
+import sys
+import os
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+
+# Import Paper's NumPy-compatible API
+from paper import numpy_api as pnp
+import numpy as np
+
+def example_basic_operations():
+    """Demonstrate basic array creation and operations."""
+    print("=" * 60)
+    print("Example 1: Basic Array Creation and Operations")
+    print("=" * 60)
+    
+    # Create arrays from data (similar to NumPy)
+    a = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)
+    b = pnp.array([[7, 8, 9], [10, 11, 12]], dtype=np.float32)
+    
+    print(f"Array a: shape={a.shape}, dtype={a.dtype}")
+    print(f"Array b: shape={b.shape}, dtype={b.dtype}")
+    
+    # Element-wise addition (lazy evaluation)
+    c = a + b
+    print(f"\nLazy addition c = a + b: {c}")
+    print(f"Is lazy? {c._is_lazy}")
+    
+    # Compute the result
+    result = c.compute()
+    print(f"\nComputed result:")
+    print(result.to_numpy())
+    
+    print()
+
+def example_scalar_multiplication():
+    """Demonstrate scalar multiplication."""
+    print("=" * 60)
+    print("Example 2: Scalar Multiplication")
+    print("=" * 60)
+    
+    a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+    
+    # Scalar multiplication (works both ways)
+    c1 = a * 2
+    c2 = 2 * a
+    
+    print("Array a:")
+    print(a.to_numpy())
+    
+    print("\nResult of a * 2:")
+    print(c1.compute().to_numpy())
+    
+    print("\nResult of 2 * a:")
+    print(c2.compute().to_numpy())
+    
+    print()
+
+def example_matrix_multiplication():
+    """Demonstrate matrix multiplication."""
+    print("=" * 60)
+    print("Example 3: Matrix Multiplication")
+    print("=" * 60)
+    
+    # Create matrices for multiplication
+    a = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)  # 2x3
+    b = pnp.array([[7, 8], [9, 10], [11, 12]], dtype=np.float32)  # 3x2
+    
+    print(f"Matrix a: shape={a.shape}")
+    print(a.to_numpy())
+    
+    print(f"\nMatrix b: shape={b.shape}")
+    print(b.to_numpy())
+    
+    # Matrix multiplication using @ operator
+    c = a @ b
+    print(f"\nResult of a @ b: shape={c.shape}")
+    print(c.compute().to_numpy())
+    
+    print()
+
+def example_chained_operations():
+    """Demonstrate chaining operations with lazy evaluation."""
+    print("=" * 60)
+    print("Example 4: Chained Operations (Lazy Evaluation)")
+    print("=" * 60)
+    
+    a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+    b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+    
+    # Complex expression: (a + b) * 2
+    # This creates a computation plan without executing
+    c = (a + b) * 2
+    
+    print("Expression: (a + b) * 2")
+    print(f"Is lazy? {c._is_lazy}")
+    print(f"Result shape: {c.shape}")
+    
+    # Execute the entire computation plan
+    result = c.compute()
+    print("\nComputed result:")
+    print(result.to_numpy())
+    
+    # Note: The framework automatically optimizes this!
+    # The fused kernel performs addition and scalar multiplication in one pass
+    
+    print()
+
+def example_array_creation_functions():
+    """Demonstrate various array creation functions."""
+    print("=" * 60)
+    print("Example 5: Array Creation Functions")
+    print("=" * 60)
+    
+    # Zeros array
+    zeros = pnp.zeros((3, 4))
+    print("Zeros array (3x4):")
+    print(zeros.to_numpy())
+    
+    # Ones array
+    ones = pnp.ones((2, 3))
+    print("\nOnes array (2x3):")
+    print(ones.to_numpy())
+    
+    # Identity matrix
+    identity = pnp.eye(4)
+    print("\nIdentity matrix (4x4):")
+    print(identity.to_numpy())
+    
+    # Random array
+    random = pnp.random_rand((2, 2))
+    print("\nRandom array (2x2) - values in [0, 1):")
+    print(random.to_numpy())
+    
+    print()
+
+def example_transpose():
+    """Demonstrate transpose operation."""
+    print("=" * 60)
+    print("Example 6: Transpose")
+    print("=" * 60)
+    
+    a = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)
+    
+    print("Original array (2x3):")
+    print(a.to_numpy())
+    
+    # Transpose using .T property
+    a_t = a.T
+    
+    print(f"\nTransposed array (3x2):")
+    print(a_t.to_numpy())
+    
+    print()
+
+def example_file_operations():
+    """Demonstrate saving and loading arrays."""
+    print("=" * 60)
+    print("Example 7: File I/O Operations")
+    print("=" * 60)
+    
+    # Create temporary directory for demo
+    import tempfile
+    temp_dir = tempfile.mkdtemp()
+    
+    # Create and save an array
+    a = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)
+    save_path = os.path.join(temp_dir, "my_array.bin")
+    
+    print(f"Saving array to: {save_path}")
+    pnp.save(save_path, a)
+    print("Array saved successfully!")
+    
+    # Load the array
+    print(f"\nLoading array from: {save_path}")
+    loaded = pnp.load(save_path, shape=(2, 3))
+    print("Loaded array:")
+    print(loaded.to_numpy())
+    
+    # Cleanup
+    import shutil
+    shutil.rmtree(temp_dir)
+    print("\nCleanup completed.")
+    
+    print()
+
+def example_large_arrays():
+    """Demonstrate out-of-core capabilities with larger arrays."""
+    print("=" * 60)
+    print("Example 8: Large Arrays (Out-of-Core)")
+    print("=" * 60)
+    
+    print("Creating large random arrays (1000x1000)...")
+    a = pnp.random_rand((1000, 1000))
+    b = pnp.random_rand((1000, 1000))
+    
+    print(f"Array a: shape={a.shape}, size={a.size} elements")
+    print(f"Array b: shape={b.shape}, size={b.size} elements")
+    
+    # Perform operations without loading entire arrays into memory
+    print("\nPerforming lazy operations...")
+    c = (a + b) * 0.5
+    
+    print(f"Result shape: {c.shape}")
+    print("Note: The operations are lazy - no computation has happened yet!")
+    
+    print("\nComputing result (this actually executes the operations)...")
+    result = c.compute()
+    print(f"Result computed successfully! Shape: {result.shape}")
+    print("The framework handled this efficiently using disk-backed matrices.")
+    
+    print()
+
+def example_numpy_comparison():
+    """Compare Paper API with NumPy syntax."""
+    print("=" * 60)
+    print("Example 9: NumPy Compatibility Comparison")
+    print("=" * 60)
+    
+    print("NumPy syntax:")
+    print("  import numpy as np")
+    print("  a = np.array([[1, 2], [3, 4]])")
+    print("  b = np.array([[5, 6], [7, 8]])")
+    print("  c = (a + b) * 2")
+    print()
+    
+    print("Paper's NumPy-compatible API syntax:")
+    print("  from paper import numpy_api as pnp")
+    print("  a = pnp.array([[1, 2], [3, 4]])")
+    print("  b = pnp.array([[5, 6], [7, 8]])")
+    print("  c = (a + b) * 2")
+    print("  result = c.compute()  # Execute the lazy computation")
+    print()
+    
+    print("The API is nearly identical, making migration easy!")
+    print("Benefits:")
+    print("  - Familiar NumPy interface")
+    print("  - Handles datasets larger than memory")
+    print("  - Automatic optimization (fusion, caching)")
+    print("  - Lazy evaluation for efficiency")
+    
+    print()
+
+def main():
+    """Run all examples."""
+    print("\n" + "=" * 60)
+    print("Paper Framework - NumPy-Compatible API Examples")
+    print("=" * 60 + "\n")
+    
+    example_basic_operations()
+    example_scalar_multiplication()
+    example_matrix_multiplication()
+    example_chained_operations()
+    example_array_creation_functions()
+    example_transpose()
+    example_file_operations()
+    example_large_arrays()
+    example_numpy_comparison()
+    
+    print("=" * 60)
+    print("All examples completed successfully!")
+    print("=" * 60)
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/simple_demo.py b/examples/simple_demo.py
new file mode 100644
index 0000000..32b4529
--- /dev/null
+++ b/examples/simple_demo.py
@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+"""
+Simple NumPy API Demo
+
+A minimal example showing how to use Paper's NumPy-compatible API.
+"""
+
+import sys
+import os
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+
+from paper import numpy_api as pnp
+import numpy as np
+
+# Create arrays using familiar NumPy syntax
+a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+
+print("Array a:")
+print(a.to_numpy())
+
+print("\nArray b:")
+print(b.to_numpy())
+
+# Build a computation plan (lazy evaluation)
+result = (a + b) * 2
+
+print(f"\nComputation plan built: {result}")
+print(f"Shape: {result.shape}")
+
+# Execute the plan
+computed = result.compute()
+
+print("\nComputed result:")
+print(computed.to_numpy())
+
+print("\n✓ Success! The NumPy-compatible API works seamlessly.")
diff --git a/paper/__init__.py b/paper/__init__.py
index 245eb99..6efa7c6 100644
--- a/paper/__init__.py
+++ b/paper/__init__.py
@@ -1,2 +1,5 @@
 from .core import PaperMatrix
-from .plan import Plan, EagerNode
\ No newline at end of file
+from .plan import Plan, EagerNode
+
+# NumPy-compatible API
+from . import numpy_api as np
\ No newline at end of file
diff --git a/paper/numpy_api.py b/paper/numpy_api.py
new file mode 100644
index 0000000..dfc83a5
--- /dev/null
+++ b/paper/numpy_api.py
@@ -0,0 +1,402 @@
+"""
+NumPy-compatible API Layer for Paper Framework
+
+This module provides a NumPy-compatible interface for the Paper framework,
+allowing users to leverage familiar NumPy operations while benefiting from
+out-of-core matrix operations for datasets larger than memory.
+
+The API mimics NumPy's interface but operates on disk-backed matrices through
+the Paper framework's lazy evaluation and optimization capabilities.
+"""
+
+import numpy as np
+import os
+import tempfile
+from typing import Union, Tuple, Optional
+from .core import PaperMatrix
+from .plan import Plan, EagerNode
+from .config import TILE_SIZE
+
+
+class ndarray:
+    """
+    NumPy-compatible array class for out-of-core matrix operations.
+    
+    This class provides a familiar NumPy-like interface while delegating
+    operations to the Paper framework's lazy evaluation system.
+    
+    Attributes:
+        shape: Tuple representing the dimensions of the array
+        dtype: NumPy data type of the array elements
+        ndim: Number of dimensions (always 2 for matrices)
+        size: Total number of elements in the array
+    """
+    
+    def __init__(self, data=None, filepath=None, shape=None, dtype=np.float32, mode='r'):
+        """
+        Initialize a Paper ndarray.
+        
+        Args:
+            data: Initial data (numpy array or nested lists)
+            filepath: Path to existing matrix file on disk
+            shape: Shape tuple for the array
+            dtype: Data type for array elements
+            mode: File access mode ('r', 'w+', 'r+')
+        """
+        self.dtype = np.dtype(dtype)
+        
+        if filepath is not None:
+            # Load from existing file
+            if shape is None:
+                raise ValueError("Shape must be provided when loading from filepath")
+            self._matrix = PaperMatrix(filepath, shape, dtype=self.dtype, mode=mode)
+            self._plan = Plan(EagerNode(self._matrix))
+            self._filepath = filepath
+            self._is_lazy = False
+        elif data is not None:
+            # Create from data
+            data_array = np.asarray(data, dtype=self.dtype)
+            if data_array.ndim != 2:
+                raise ValueError("Paper arrays must be 2-dimensional matrices")
+            
+            # Create a temporary file to store the data
+            fd, self._filepath = tempfile.mkstemp(suffix='.bin', prefix='paper_array_')
+            os.close(fd)
+            
+            # Write data to file
+            self._matrix = PaperMatrix(self._filepath, data_array.shape, dtype=self.dtype, mode='w+')
+            self._matrix.data[:] = data_array
+            self._matrix.data.flush()
+            self._plan = Plan(EagerNode(self._matrix))
+            self._is_lazy = False
+        elif shape is not None:
+            # Create empty array
+            fd, self._filepath = tempfile.mkstemp(suffix='.bin', prefix='paper_array_')
+            os.close(fd)
+            
+            self._matrix = PaperMatrix(self._filepath, shape, dtype=self.dtype, mode='w+')
+            self._plan = Plan(EagerNode(self._matrix))
+            self._is_lazy = False
+        else:
+            raise ValueError("Must provide either data, filepath, or shape")
+    
+    @classmethod
+    def _from_plan(cls, plan: Plan, shape: Tuple[int, int], dtype=np.float32):
+        """Internal constructor for lazy arrays from computation plans."""
+        obj = cls.__new__(cls)
+        obj._plan = plan
+        obj._matrix = None
+        obj._filepath = None
+        obj.dtype = np.dtype(dtype)
+        obj._is_lazy = True
+        obj._shape = shape
+        return obj
+    
+    @property
+    def shape(self) -> Tuple[int, int]:
+        """Return the shape of the array."""
+        if self._is_lazy:
+            return self._shape
+        return self._matrix.shape if self._matrix else self._plan.shape
+    
+    @property
+    def ndim(self) -> int:
+        """Return the number of dimensions (always 2 for matrices)."""
+        return 2
+    
+    @property
+    def size(self) -> int:
+        """Return the total number of elements."""
+        return self.shape[0] * self.shape[1]
+    
+    @property
+    def T(self) -> 'ndarray':
+        """
+        Return the transpose of the array.
+        Note: Currently creates a copy. Future optimization: lazy transpose.
+        """
+        # For now, materialize and transpose
+        result = self._materialize()
+        transposed = result.T
+        return array(transposed, dtype=self.dtype)
+    
+    def _materialize(self) -> np.ndarray:
+        """
+        Materialize the lazy computation into an actual NumPy array.
+        Warning: This loads the entire array into memory.
+        
+        Internal method. Use to_numpy() for public API.
+        """
+        if self._is_lazy:
+            # Compute the plan
+            fd, temp_path = tempfile.mkstemp(suffix='.bin', prefix='paper_materialized_')
+            os.close(fd)
+            
+            result_matrix, _ = self._plan.compute(temp_path)
+            data = np.array(result_matrix.data, copy=True)
+            result_matrix.close()
+            os.unlink(temp_path)
+            return data
+        else:
+            # Already materialized
+            return np.array(self._matrix.data, copy=True)
+    
+    def to_numpy(self) -> np.ndarray:
+        """
+        Convert the array to a NumPy array.
+        
+        This method materializes the array, loading it into memory.
+        For lazy arrays, this executes the computation plan first.
+        
+        Returns:
+            np.ndarray: A NumPy array containing the data
+            
+        Warning:
+            This loads the entire array into memory and may fail for
+            very large arrays that don't fit in RAM.
+        
+        Examples:
+            >>> import paper.numpy_api as pnp
+            >>> a = pnp.array([[1, 2], [3, 4]])
+            >>> numpy_arr = a.to_numpy()
+            >>> print(numpy_arr)
+            [[1. 2.]
+             [3. 4.]]
+        """
+        return self._materialize()
+    
+    def compute(self, output_path: Optional[str] = None, cache_size_tiles: Optional[int] = None):
+        """
+        Execute the lazy computation plan and return a materialized ndarray.
+        
+        Args:
+            output_path: Optional path to save the result
+            cache_size_tiles: Optional cache size for buffer manager
+            
+        Returns:
+            ndarray: Materialized result
+        """
+        if not self._is_lazy:
+            return self
+        
+        if output_path is None:
+            fd, output_path = tempfile.mkstemp(suffix='.bin', prefix='paper_computed_')
+            os.close(fd)
+        
+        result_matrix, _ = self._plan.compute(output_path, cache_size_tiles)
+        
+        # Create a new ndarray from the result
+        result = ndarray.__new__(ndarray)
+        result._matrix = result_matrix
+        result._filepath = output_path
+        result._plan = Plan(EagerNode(result_matrix))
+        result.dtype = self.dtype
+        result._is_lazy = False
+        return result
+    
+    def __add__(self, other: Union['ndarray', int, float]) -> 'ndarray':
+        """Element-wise addition."""
+        if isinstance(other, (int, float)):
+            # Scalar addition - not yet optimized
+            raise NotImplementedError("Scalar addition not yet implemented")
+        
+        if not isinstance(other, ndarray):
+            raise TypeError(f"Unsupported operand type for +: 'ndarray' and '{type(other).__name__}'")
+        
+        if self.shape != other.shape:
+            raise ValueError(f"Shape mismatch: {self.shape} vs {other.shape}")
+        
+        # Create lazy addition plan
+        new_plan = self._plan + other._plan
+        return ndarray._from_plan(new_plan, self.shape, self.dtype)
+    
+    def __mul__(self, other: Union['ndarray', int, float]) -> 'ndarray':
+        """Element-wise multiplication or scalar multiplication."""
+        if isinstance(other, (int, float)):
+            # Scalar multiplication
+            new_plan = self._plan * other
+            return ndarray._from_plan(new_plan, self.shape, self.dtype)
+        
+        # Element-wise multiplication not yet implemented
+        raise NotImplementedError("Element-wise multiplication with arrays not yet implemented")
+    
+    def __rmul__(self, other: Union[int, float]) -> 'ndarray':
+        """Right scalar multiplication."""
+        return self.__mul__(other)
+    
+    def __matmul__(self, other: 'ndarray') -> 'ndarray':
+        """Matrix multiplication."""
+        if not isinstance(other, ndarray):
+            raise TypeError(f"Unsupported operand type for @: 'ndarray' and '{type(other).__name__}'")
+        
+        if self.shape[1] != other.shape[0]:
+            raise ValueError(f"Inner dimensions must match: {self.shape} @ {other.shape}")
+        
+        # Create lazy matmul plan
+        new_plan = self._plan @ other._plan
+        result_shape = (self.shape[0], other.shape[1])
+        return ndarray._from_plan(new_plan, result_shape, self.dtype)
+    
+    def __repr__(self) -> str:
+        """String representation."""
+        if self._is_lazy:
+            return f"ndarray(shape={self.shape}, dtype={self.dtype.name}, lazy=True)"
+        return f"ndarray(shape={self.shape}, dtype={self.dtype.name})"
+    
+    def __del__(self):
+        """Cleanup temporary files."""
+        if hasattr(self, '_matrix') and self._matrix is not None:
+            try:
+                self._matrix.close()
+            except:
+                pass
+        
+        # Clean up temporary files
+        if hasattr(self, '_filepath') and self._filepath and os.path.exists(self._filepath):
+            if 'paper_array_' in self._filepath or 'paper_computed_' in self._filepath:
+                try:
+                    os.unlink(self._filepath)
+                except:
+                    pass
+
+
+def array(data, dtype=np.float32) -> ndarray:
+    """
+    Create a Paper array from existing data.
+    
+    Args:
+        data: Array-like data (list, tuple, or numpy array)
+        dtype: Data type for the array elements
+        
+    Returns:
+        ndarray: A Paper ndarray object
+        
+    Examples:
+        >>> import paper.numpy_api as pnp
+        >>> a = pnp.array([[1, 2], [3, 4]])
+        >>> b = pnp.array([[5, 6], [7, 8]])
+        >>> c = a + b  # Lazy operation
+        >>> result = c.compute()  # Execute
+    """
+    return ndarray(data=data, dtype=dtype)
+
+
+def zeros(shape: Tuple[int, int], dtype=np.float32) -> ndarray:
+    """
+    Create an array filled with zeros.
+    
+    Args:
+        shape: Shape of the array (rows, cols)
+        dtype: Data type
+        
+    Returns:
+        ndarray: Array filled with zeros
+    """
+    arr = ndarray(shape=shape, dtype=dtype)
+    arr._matrix.data[:] = 0
+    arr._matrix.data.flush()
+    return arr
+
+
+def ones(shape: Tuple[int, int], dtype=np.float32) -> ndarray:
+    """
+    Create an array filled with ones.
+    
+    Args:
+        shape: Shape of the array (rows, cols)
+        dtype: Data type
+        
+    Returns:
+        ndarray: Array filled with ones
+    """
+    arr = ndarray(shape=shape, dtype=dtype)
+    arr._matrix.data[:] = 1
+    arr._matrix.data.flush()
+    return arr
+
+
+def random_rand(shape: Tuple[int, int], dtype=np.float32) -> ndarray:
+    """
+    Create an array with random values from [0, 1).
+    
+    Args:
+        shape: Shape of the array (rows, cols)
+        dtype: Data type
+        
+    Returns:
+        ndarray: Array with random values
+    """
+    arr = ndarray(shape=shape, dtype=dtype)
+    
+    # Fill tile by tile to avoid loading entire matrix into memory
+    for r_start in range(0, shape[0], TILE_SIZE):
+        r_end = min(r_start + TILE_SIZE, shape[0])
+        for c_start in range(0, shape[1], TILE_SIZE):
+            c_end = min(c_start + TILE_SIZE, shape[1])
+            tile_shape = (r_end - r_start, c_end - c_start)
+            random_tile = np.random.rand(*tile_shape).astype(dtype)
+            arr._matrix.data[r_start:r_end, c_start:c_end] = random_tile
+    
+    arr._matrix.data.flush()
+    return arr
+
+
+def eye(n: int, dtype=np.float32) -> ndarray:
+    """
+    Create a 2-D identity matrix.
+    
+    Args:
+        n: Number of rows and columns
+        dtype: Data type
+        
+    Returns:
+        ndarray: Identity matrix
+    """
+    arr = zeros((n, n), dtype=dtype)
+    
+    # Set diagonal to 1
+    for i in range(n):
+        arr._matrix.data[i, i] = 1
+    
+    arr._matrix.data.flush()
+    return arr
+
+
+def load(filepath: str, shape: Tuple[int, int], dtype=np.float32) -> ndarray:
+    """
+    Load a matrix from a binary file.
+    
+    Args:
+        filepath: Path to the binary matrix file
+        shape: Shape of the matrix
+        dtype: Data type
+        
+    Returns:
+        ndarray: Loaded array
+    """
+    return ndarray(filepath=filepath, shape=shape, dtype=dtype, mode='r')
+
+
+def save(filepath: str, arr: ndarray):
+    """
+    Save an array to a binary file.
+    
+    Args:
+        filepath: Path to save the file
+        arr: Array to save
+    """
+    if arr._is_lazy:
+        # Compute first if lazy
+        arr = arr.compute(output_path=filepath)
+    else:
+        # Copy to new location if needed
+        if arr._filepath != filepath:
+            # Simple file copy
+            import shutil
+            shutil.copy(arr._filepath, filepath)
+
+
+# Expose common NumPy functions
+dot = lambda a, b: a @ b  # Matrix multiplication
+add = lambda a, b: a + b  # Element-wise addition
+multiply = lambda a, b: a * b  # Scalar or element-wise multiplication
diff --git a/paper/plan.py b/paper/plan.py
index 223242f..47599ea 100644
--- a/paper/plan.py
+++ b/paper/plan.py
@@ -161,6 +161,8 @@ def execute(self, output_path, buffer_manager):
                     C.data[r:r_end, c:c_end] = TMP.data[r:r_end, c:c_end] * self.right
             
             C.data.flush()
-            TMP.close()
+            # Only close TMP if it's not from an EagerNode (i.e., it was computed, not an input)
+            if not isinstance(self.left, EagerNode):
+                TMP.close()
             return C
 
diff --git a/tests/test_numpy_api.py b/tests/test_numpy_api.py
new file mode 100644
index 0000000..315f89b
--- /dev/null
+++ b/tests/test_numpy_api.py
@@ -0,0 +1,417 @@
+"""
+Unit tests for the NumPy-compatible API layer.
+
+Tests the numpy_api module to ensure it provides a familiar NumPy interface
+while leveraging the Paper framework's out-of-core capabilities.
+"""
+
+import unittest
+import os
+import tempfile
+import shutil
+import numpy as np
+import sys
+
+# Add the parent directory to the path so we can import the paper module
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+
+from paper import numpy_api as pnp
+
+
+class TestNumpyAPIArrayCreation(unittest.TestCase):
+    """Test cases for array creation functions."""
+    
+    def setUp(self):
+        """Set up test fixtures."""
+        self.test_dir = tempfile.mkdtemp()
+    
+    def tearDown(self):
+        """Clean up after tests."""
+        if os.path.exists(self.test_dir):
+            shutil.rmtree(self.test_dir)
+    
+    def test_array_creation_from_list(self):
+        """Test creating array from nested lists."""
+        data = [[1, 2, 3], [4, 5, 6]]
+        arr = pnp.array(data)
+        
+        self.assertEqual(arr.shape, (2, 3))
+        self.assertEqual(arr.dtype, np.float32)
+        self.assertEqual(arr.ndim, 2)
+        self.assertEqual(arr.size, 6)
+    
+    def test_array_creation_from_numpy_array(self):
+        """Test creating array from numpy array."""
+        np_arr = np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float32)
+        arr = pnp.array(np_arr)
+        
+        self.assertEqual(arr.shape, (2, 2))
+        self.assertEqual(arr.dtype, np.float32)
+        
+        # Verify data correctness
+        materialized = arr.to_numpy()
+        np.testing.assert_array_equal(materialized, np_arr)
+    
+    def test_array_creation_with_dtype(self):
+        """Test array creation with different dtypes."""
+        for dtype in [np.float32, np.float64]:
+            with self.subTest(dtype=dtype):
+                arr = pnp.array([[1, 2], [3, 4]], dtype=dtype)
+                self.assertEqual(arr.dtype, dtype)
+    
+    def test_zeros_creation(self):
+        """Test zeros array creation."""
+        arr = pnp.zeros((3, 4))
+        
+        self.assertEqual(arr.shape, (3, 4))
+        self.assertEqual(arr.dtype, np.float32)
+        
+        # Verify all zeros
+        materialized = arr.to_numpy()
+        np.testing.assert_array_equal(materialized, np.zeros((3, 4), dtype=np.float32))
+    
+    def test_ones_creation(self):
+        """Test ones array creation."""
+        arr = pnp.ones((2, 5))
+        
+        self.assertEqual(arr.shape, (2, 5))
+        
+        # Verify all ones
+        materialized = arr.to_numpy()
+        np.testing.assert_array_equal(materialized, np.ones((2, 5), dtype=np.float32))
+    
+    def test_random_rand_creation(self):
+        """Test random array creation."""
+        arr = pnp.random_rand((3, 3))
+        
+        self.assertEqual(arr.shape, (3, 3))
+        
+        # Verify values are in [0, 1)
+        materialized = arr.to_numpy()
+        self.assertTrue(np.all(materialized >= 0))
+        self.assertTrue(np.all(materialized < 1))
+    
+    def test_eye_creation(self):
+        """Test identity matrix creation."""
+        arr = pnp.eye(4)
+        
+        self.assertEqual(arr.shape, (4, 4))
+        
+        # Verify identity matrix
+        materialized = arr.to_numpy()
+        np.testing.assert_array_equal(materialized, np.eye(4, dtype=np.float32))
+    
+    def test_array_properties(self):
+        """Test array properties (shape, ndim, size, dtype)."""
+        arr = pnp.array([[1, 2, 3], [4, 5, 6]])
+        
+        self.assertEqual(arr.shape, (2, 3))
+        self.assertEqual(arr.ndim, 2)
+        self.assertEqual(arr.size, 6)
+        self.assertEqual(arr.dtype, np.float32)
+    
+    def test_to_numpy_method(self):
+        """Test to_numpy() method for converting to NumPy arrays."""
+        # Test with materialized array
+        arr = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        numpy_arr = arr.to_numpy()
+        
+        self.assertIsInstance(numpy_arr, np.ndarray)
+        expected = np.array([[1, 2], [3, 4]], dtype=np.float32)
+        np.testing.assert_array_equal(numpy_arr, expected)
+        
+        # Test with lazy array
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+        c = a + b
+        
+        self.assertTrue(c._is_lazy)
+        numpy_result = c.to_numpy()
+        
+        self.assertIsInstance(numpy_result, np.ndarray)
+        expected_lazy = np.array([[6, 8], [10, 12]], dtype=np.float32)
+        np.testing.assert_array_equal(numpy_result, expected_lazy)
+
+
+class TestNumpyAPIOperations(unittest.TestCase):
+    """Test cases for NumPy-compatible operations."""
+    
+    def setUp(self):
+        """Set up test fixtures."""
+        self.test_dir = tempfile.mkdtemp()
+    
+    def tearDown(self):
+        """Clean up after tests."""
+        if os.path.exists(self.test_dir):
+            shutil.rmtree(self.test_dir)
+    
+    def test_addition_operation(self):
+        """Test array addition."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+        
+        c = a + b
+        
+        # Check that result is lazy
+        self.assertTrue(c._is_lazy)
+        self.assertEqual(c.shape, (2, 2))
+        
+        # Compute and verify
+        result = c.compute()
+        expected = np.array([[6, 8], [10, 12]], dtype=np.float32)
+        
+        materialized = result.to_numpy()
+        np.testing.assert_array_equal(materialized, expected)
+    
+    def test_scalar_multiplication(self):
+        """Test scalar multiplication."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        
+        # Test both left and right multiplication
+        c1 = a * 2
+        c2 = 2 * a
+        
+        # Both should be lazy
+        self.assertTrue(c1._is_lazy)
+        self.assertTrue(c2._is_lazy)
+        
+        # Compute and verify
+        result1 = c1.compute()
+        result2 = c2.compute()
+        expected = np.array([[2, 4], [6, 8]], dtype=np.float32)
+        
+        np.testing.assert_array_equal(result1.to_numpy(), expected)
+        np.testing.assert_array_equal(result2.to_numpy(), expected)
+    
+    def test_matrix_multiplication(self):
+        """Test matrix multiplication."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+        
+        c = a @ b
+        
+        # Check lazy and shape
+        self.assertTrue(c._is_lazy)
+        self.assertEqual(c.shape, (2, 2))
+        
+        # Compute and verify
+        result = c.compute()
+        expected = np.array([[1, 2], [3, 4]], dtype=np.float32) @ np.array([[5, 6], [7, 8]], dtype=np.float32)
+        
+        materialized = result.to_numpy()
+        np.testing.assert_array_almost_equal(materialized, expected, decimal=5)
+    
+    def test_chained_operations(self):
+        """Test chaining multiple operations."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+        
+        # (a + b) * 2
+        c = (a + b) * 2
+        
+        self.assertTrue(c._is_lazy)
+        
+        # Compute and verify
+        result = c.compute()
+        expected = (np.array([[1, 2], [3, 4]], dtype=np.float32) + 
+                   np.array([[5, 6], [7, 8]], dtype=np.float32)) * 2
+        
+        materialized = result.to_numpy()
+        np.testing.assert_array_equal(materialized, expected)
+    
+    def test_matrix_multiplication_different_sizes(self):
+        """Test matrix multiplication with non-square matrices."""
+        a = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)  # 2x3
+        b = pnp.array([[7, 8], [9, 10], [11, 12]], dtype=np.float32)  # 3x2
+        
+        c = a @ b
+        
+        self.assertEqual(c.shape, (2, 2))
+        
+        # Compute and verify
+        result = c.compute()
+        expected = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32) @ np.array([[7, 8], [9, 10], [11, 12]], dtype=np.float32)
+        
+        materialized = result.to_numpy()
+        np.testing.assert_array_almost_equal(materialized, expected, decimal=5)
+    
+    def test_transpose_operation(self):
+        """Test array transpose."""
+        a = pnp.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)
+        
+        b = a.T
+        
+        self.assertEqual(b.shape, (3, 2))
+        
+        # Verify transpose correctness
+        materialized = b.to_numpy()
+        expected = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32).T
+        
+        np.testing.assert_array_equal(materialized, expected)
+
+
+class TestNumpyAPIErrorHandling(unittest.TestCase):
+    """Test cases for error handling in NumPy API."""
+    
+    def test_shape_mismatch_addition(self):
+        """Test that addition with mismatched shapes raises error."""
+        a = pnp.array([[1, 2], [3, 4]])
+        b = pnp.array([[1, 2, 3]])
+        
+        with self.assertRaises(ValueError):
+            c = a + b
+    
+    def test_dimension_mismatch_matmul(self):
+        """Test that matmul with incompatible dimensions raises error."""
+        a = pnp.array([[1, 2], [3, 4]])  # 2x2
+        b = pnp.array([[1, 2, 3]])  # 1x3
+        
+        with self.assertRaises(ValueError):
+            c = a @ b
+    
+    def test_invalid_array_creation(self):
+        """Test that creating array without data or shape raises error."""
+        with self.assertRaises(ValueError):
+            arr = pnp.ndarray()
+
+
+class TestNumpyAPIFileOperations(unittest.TestCase):
+    """Test cases for file I/O operations."""
+    
+    def setUp(self):
+        """Set up test fixtures."""
+        self.test_dir = tempfile.mkdtemp()
+    
+    def tearDown(self):
+        """Clean up after tests."""
+        if os.path.exists(self.test_dir):
+            shutil.rmtree(self.test_dir)
+    
+    def test_load_array(self):
+        """Test loading array from file."""
+        # Create a test file
+        test_path = os.path.join(self.test_dir, "test.bin")
+        test_data = np.array([[1, 2], [3, 4]], dtype=np.float32)
+        test_data.tofile(test_path)
+        
+        # Load using Paper API
+        arr = pnp.load(test_path, shape=(2, 2))
+        
+        self.assertEqual(arr.shape, (2, 2))
+        
+        # Verify data
+        materialized = arr.to_numpy()
+        np.testing.assert_array_equal(materialized, test_data)
+    
+    def test_save_array(self):
+        """Test saving array to file."""
+        arr = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        
+        save_path = os.path.join(self.test_dir, "saved.bin")
+        pnp.save(save_path, arr)
+        
+        self.assertTrue(os.path.exists(save_path))
+        
+        # Load and verify
+        loaded_data = np.fromfile(save_path, dtype=np.float32).reshape((2, 2))
+        expected = np.array([[1, 2], [3, 4]], dtype=np.float32)
+        np.testing.assert_array_equal(loaded_data, expected)
+    
+    def test_save_lazy_array(self):
+        """Test saving lazy (uncomputed) array."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+        
+        c = a + b  # Lazy
+        
+        save_path = os.path.join(self.test_dir, "lazy_saved.bin")
+        pnp.save(save_path, c)
+        
+        self.assertTrue(os.path.exists(save_path))
+        
+        # Load and verify
+        loaded_data = np.fromfile(save_path, dtype=np.float32).reshape((2, 2))
+        expected = np.array([[6, 8], [10, 12]], dtype=np.float32)
+        np.testing.assert_array_equal(loaded_data, expected)
+
+
+class TestNumpyAPIHelperFunctions(unittest.TestCase):
+    """Test cases for helper functions."""
+    
+    def test_dot_function(self):
+        """Test dot product (matrix multiplication) function."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+        
+        c = pnp.dot(a, b)
+        
+        # Should be same as a @ b
+        self.assertTrue(c._is_lazy)
+        self.assertEqual(c.shape, (2, 2))
+    
+    def test_add_function(self):
+        """Test add function."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        b = pnp.array([[5, 6], [7, 8]], dtype=np.float32)
+        
+        c = pnp.add(a, b)
+        
+        # Should be same as a + b
+        self.assertTrue(c._is_lazy)
+        self.assertEqual(c.shape, (2, 2))
+    
+    def test_multiply_function(self):
+        """Test multiply function."""
+        a = pnp.array([[1, 2], [3, 4]], dtype=np.float32)
+        
+        c = pnp.multiply(a, 2)
+        
+        # Should be same as a * 2
+        self.assertTrue(c._is_lazy)
+        self.assertEqual(c.shape, (2, 2))
+
+
+class TestNumpyAPILargeArrays(unittest.TestCase):
+    """Test cases for larger arrays to verify out-of-core capabilities."""
+    
+    def setUp(self):
+        """Set up test fixtures."""
+        self.test_dir = tempfile.mkdtemp()
+    
+    def tearDown(self):
+        """Clean up after tests."""
+        if os.path.exists(self.test_dir):
+            shutil.rmtree(self.test_dir)
+    
+    def test_large_array_creation(self):
+        """Test creating and operating on larger arrays."""
+        # Create moderately large arrays (1000x1000)
+        a = pnp.random_rand((1000, 1000))
+        b = pnp.random_rand((1000, 1000))
+        
+        self.assertEqual(a.shape, (1000, 1000))
+        self.assertEqual(b.shape, (1000, 1000))
+        
+        # Test addition
+        c = a + b
+        self.assertTrue(c._is_lazy)
+        self.assertEqual(c.shape, (1000, 1000))
+    
+    def test_large_matrix_multiplication(self):
+        """Test matrix multiplication on larger matrices."""
+        # Create moderately sized matrices
+        a = pnp.random_rand((500, 600))
+        b = pnp.random_rand((600, 400))
+        
+        c = a @ b
+        
+        self.assertEqual(c.shape, (500, 400))
+        
+        # Just verify it can be computed without errors
+        # (full verification would be too slow)
+        result = c.compute()
+        self.assertEqual(result.shape, (500, 400))
+
+
+if __name__ == '__main__':
+    unittest.main()